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PR£ÊTACE 


Les possibilités offertes par les ordinateurs pour la réali- 
sation des expériences do calcul ont accéléré lo processus 
de mathématisation scientifique et technique. Un nombre de 
plus en plus grand de professions exigent des connaissances 
mathématiques approfondies. Pour subvenir aux besoins de 
l'enseignement, de nombreux mémentos et monographies ont 
été édités. Le présent ouvrage appartient à cette série. 

Cet ouvrage est consacré aux questions liées à l’élabora- 
tion des méthodes de solution numérique des problèmes de 
mécanique et de physique mathématique. Actuellement, 
la théorie est en pleine essor et est loin d'être achevée. 
Néanmoins un certain nombre d'idées et de notions assez 
simples sont devenues fondamentales, elles sont exposées 
dans le présent ouvrage. Il est évident que la connaissance 
de l'alphabet des méthodes numériques ne saurait faire 
du lecteur un calculateur qualifié. Il faut à cet effet 
une étude plus profonde de la littérature et surtout une 
grande expérience personnelle de résolution des problèmes 
concrets qui deviennent bien plus efficaces lorsque les fonde- 
ments de la théorie sont clairs. 

La diversité des professions et du niveau de préparation 
du lecteur ont obligé l’auteur à renoncer à un exposé formel, 
traditionnel pour la littérature mathématique. Il est ainsi 
possible que certains mathématiciens professionnels raffinés 
trouvent cet ouvrage banal et peut-être même vulgaire, ct 
nous sommes d'accord avec eux. Mais un « armement » 
mathématique complet, tenant compte, de toutes les éven- 
tualités logiques n'est indispensable que si l’on no dispose 
que de la logique mathématique. Mais il n’en est pas ainsi 
pour la majorité des objectifs réels, ce qui permet d'écrire des 
ouvrages tels que celui-ci. D'un autre coté nous désirons 
que ce livre soit lu, nous avons donc essayé de faire dge:telle 
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sorte que Île lecteur n'ait pas de raisons pour le 
rejeter. n 

L'exposé est mené au niveau de rigueur « physique ». 
L'auteur tient compte moins de la préparation mathémati- 
que du lecteur que de son bon sens ct de son intelligence. 
En général, l'étude de telle ou telle question est menée sur 
un exemple typique simple, puis la possibilité de généralisa- 
tion des résultats obtenus à des cas plus compliqués est 
cnvisagée. 

À la fin de chaque paragraphe on trouve des problèmes 
qui sont de difficultés différentes. La solution de la majorité 
de ces problèmes exige une compréhension claire des principes 
exposés dans le texte et la faculté de les développer indé- 
pendamment. 

Les trois premiers paragraphes servent d'introduction 
et concernent les questions classiques de calcul numérique, 
à savoir les méthodes itératives de résolution des équations, 
les formules d'interpolation cet de quadrature, l'intégration 
numérique des équations différentielles ordinaires. 

Les autres paragraphes formant le contenu essentiel de 
l'ouvrage sont consacrés aux méthodes numériques de résolu- 
tion des équations aux dérivées partielles, à l'élaboration 
et à l'étudo des algorithmes de calcul correspondants. 
On y envisage les procédés concrets d'approximation et de 
stabilité des problèmes aux différences, les propriétés des 
schémas explicites et implicites, les méthodes d'élaboration 
des formules de calcul, les méthodes de résolution des 
systèmes d'équations aux différences, la possibilité de réali- 
sation des algorithmes, etc.: 

Nous proposons donc au lecteur une théorie simplifiée des 
méthodes numériques modernes de résolution des problèmes 
de mécanique, de dynamique des gaz et de physique mathé- 
matique donnant lieu à l'intégration des équations diffé- 
renticlles. 

Cet ouvrage est le fruit de longues années de travail de 
l'auteur dans le domaine des mathématiques appliquées et de 
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conférences faites à l'Ecole Physico-technique de Moscou. 
J1 peut être utile aux étudiants des facultés des sciences et 
des Grandes Ecoles pour une première connaissance des 
méthodes de calcul numérique. 

Les questions abordées dans les trois premiers para- 
graphes sont exposées dans de nombreux ouvrages et on 
peut les trouver en détail dans tous les cours de calcul 
approché. 

Pour une étude plus approfondie des méthodes do résolu- 
lion des équations aux dérivées particiles on peut recom- 
mander les ouvrages suivants: 

1. lonyuon C.K.,PaGenLrKkniü B. C. Bsene- 
HIIC B TEOI HIO Pa8HOCTHBX CXeM (Introduction à la théorie 
des schémas aux différences). Fsmarrus, 1962. 

2. Camapckni À. À. Bgeneuxne B Teopnio pas- 
HOCTHEIX CXeM (Introduction à la théorie des schémas aux 
différences). « Hayxa », 1971. 

3. fnenko H. H. Merox npoôusx maros pewcnna 
MHOTOMEOPHBIX 8AAA4 MATeMaTudeCKOÏ busnkn (Méthode des 
pas fractionnaires de résolution des problèmes multidimen- 
sionnels de physique mathématique). « Hayxa », 1967. 

4. Richtmyeor:R. Difference Methods for Initial- 
Value Problems. New York, 1957. 

5. Wasow W., Forsythe G. Finite-difference 
Methods for Partial Differential Equations, New York, 
1960. 

L'auteur tend à remercier V. S. Riabenki pour son 
concours lors de la parution du présent ouvrage. 


V. Diatchenko 


INTRODUCTION 


On appelle méthode numérique de résolution d'un pro- 
blèmo uno certaine suite d'opérations sur des nombres, 
c'est-à-dire un algorithme de calcul dont le langage est 
formé par des chiffres ct des opérations arithmétiques. 
Ce langage primitif permet de réaliser les méthodes numéri- 
ques sur ordinateur, ce qui fait qu'elles sont un instrument 
puissant et universel de recherche. 

Cependant les problèmes à résoudre sc formulent on géné- 
ral dans le langage mathématique ordinaire (équations, 
fonctions, opérateurs différentiels, etc.). Ainsi l'élaboration 
d'une méthode numérique suppose obligatoirement la subs- 
titution, l'approximation du problème initial par un 
autre voisin formulé en termes de nombres et d'opérations 
arithmétiques. Malgré la diversité des méthodes de cette 
substitution certaines propriétés générales leur sont com- 
:munes. 

Considérons un exemple très simple. Il y} a lieu de 
trouver la solution de l'équation suivante 


x — a = 0, a > O0, (1) 


c'est-à-dire d'extraire la racine carrée d’un nombre donné 4. 
On peut évidemment écrire + = Wa, mais le symbole Y 
ne résout pas le problème, ne donne pas une méthode de 
calcul de la grandeur x. 

Nous allons procéder de la manière suivante. Nous allons 
nous donner une approximation initiale x, (par exemple 
zo = 1) et successivement à l'aide de la formule 


tn = (ani + À) (2) 


Tn-1 


calculer les valeurs z,, x, . . . Nous allons arrêter ce proces- 
sus sur un certain nr = À, lo résultat obtenu x, sera déclaré 
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être la solution approchée du problème initial (1), c'est-à- 


dire : 
V'a TN: 

La possibilité de cette hypothèse dépend évidemment des 
exigences imposées à la précision de la solution, do la gran- 
deur a ot du paramètro Ÿ. Pour des conditions quelconques, 
il y a lieu de démontrer que, pour lout a, par un choix con- 
venable de Ÿ on peut faire æ, aussi voisin que l'on veut de 
la valeur exacte Va. 

Nous allons démontrer que notre alogrithmo (2) satisfait 
à cotte condition. Re 


= = 1<Le,. (3) 


Divisions l'égalité (2) par Va et substituons (3), il vient 


1 1 
tes (1+enitr—) , 
d'où 
1 1 ef _s 


1 
past )sra.  @ 


Comme 1 + e, = 1/V a > 0, en vertu de la dernière égalité, 
tous les e, à partir du premier sont positifs. Ceci signifie que 


_En-1 
En-i +1 É 


Compte tenu de ce qui vient d'être dit, on obtient à partir 
de (4) 


Î 
En < 9 En-15 (5) 


c'est-à-dire que €; décroît lorsque n augmente, et ceci plus 
rapidement qu'une progression géométrique de raison 1/2. 
Par conséquent 


z, Va pour V —+ oo, (6) 
notre assertion so trouve ainsi démontrée. 
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Fig. 1. 


Sur la fig. 1 on peut voir le processus itératif (2). Les 
graphiques correspondant aux membres premier y, (x) et 


second y, (x) de (2). Comme de toute évidence y, (Va) — 


— y, (Va), ces graphiques se rencontrent au point x = Wa. 
Les itérations correspondant à la formule (2) équivalent 
à un mouvement le long de la ligne brisée représentée sur la 
figure entre y, (x) et y: (x). On peut ainsi se rendre compte 


une fois de plus de ce que lesitérations convergent vers Va, 
pour V — co. 

Lors de l'étude de la convergence nous avons en une cer- 
taine mesure idéalisé l'algorithme en supposant que le 
calcul exact suivant la formule (2) est possible, mais ni 
l'homme, ni l'ordinateur ne peuvent opérer avec des nombres 
réels quelconques, les calculs s'effectuent toujours avec 
un nombre restreint de décimales ct la précision du résultat 
ne saurait être supérieure à la précision du calcul. Il importe 
donc d'établir la relation existant entre ces précisions, de 
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voir si les erreurs apparaissant lors de l'arrondissement ne 
s'accumulent pas, ce qui aurait pour effet de donner un 
résultat dépourvu de toute valeur. 

Bien qu'il soit presque évident que, dans notre exemple, 
ceci n'aura pas lieu, nous allons procéder à la vérification 
formelle. L'arrondissement revient en fait à remplacer la 
formule (2) par la formule 


En — _ (an +=—) (1 na On); (7) 


a 
Tn—1 


où le facteur 1 + 6, tient compte de l'erreur d'arrondisse- 


ment introduite lors du #-ième pas du calcul, +, étant la 
suite obtenue en fait. La grandeur ô € 1 caractérise la pré- 


cision des calculs. En remplaçant x, par Va (1 + e,) on 
obtient au lieu de (4) 
2 
14 En- 
En = HU + Ôn) + ôn. 

On peut voir qu'avec l'augmentation de x e, décroît jusqu'à 
une grandeur de l'ordre de ô,, c'est-à-dire que la précision 
du résultat correspond à la précision des calculs. 

Malgré sa simplicité, l'exemple envisagé met en évidence 
d'une manière très claire les principes généraux suivants 
caractérisant toutes les méthodes de calcul. 

Primo, le probléme initial (1) est remplacé par un autre 
problème, à savoir un algorithme de calcul (2). 

Secondo, le problème (2) contient un paramètre N ne 
figurant pas dans le problème initial. 

Tertio, par un choix convenable de ce paramètre on peut 
en principe rendre la solution x, du second problème aussi 


voisine que l'on veut de la solution du premier Va. 

Enfin quarto, si l'algorithme est réalisé avec une préci- 
sion insuffisante, vu l'arrondissement, ceci n'a en fait pas 
d'influence sur ses propriétés. 


Chapitre 
premier 


CALCUL DES RACINES 
S + D'UNE ÉQUATION 


Soit à trouver la racinc réelle de l'équation 
f (x) = 0. (8) 


Nous allons supposer que cette racine se trouve à l'intérieur 
d'un certain intervalle donné (celui-ci peut être grand). 
Lo problème (1) étudié ci-dessus est un cas particulier 
de celui onvisagé ici. La méthode d'’itérations utilisée peut 
êtro étendue au cas général de l'équation (8). À cet cffet 
il y a lieu de réaliser un processus ilératif du type 
Tn —= (Tn-1) (9) 
(comparer avec (2)) convergeant vers la racine de l'équation 
(8) que nous désignerons par À. 

Nous allons voir les conditions auxquelles doit satisfaire 
la fonction œ (x) et certaines méthodes permettant de la 
construire. 

Supposons que le processus d'ilérations (9) "converge, 
c'est-à-dire quo x, —> À pour #2 —+> co. En passant à la 
limite dans les premier el second membres de (9), on obtient 
X = p(X), soil x = X doit êlre une racine commune à 
l'équation (8) et à 


x —= (x). (10) 


Donc pour obtenir la formule d'itérations, il suffit simple- 
ment d'écrire l'équation (8) sous la forme (10) (comparer (1) 
et (2)), ce que l'on peut de toute évidenco faire de différentes 
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manières. Ainsi l'équation (1) pout, en plus de la forme (2), 
s'écrire comme suit: 


a 
LL = — 
T 
ou 
a 
T=ÈÛr— —. 
A] 


II est facile de voir que ni l’une ni l’autre de ces formules 
ne peut être utilisée pour des itérations. La première donne 
di = Alto, Lo = Lo, Ta = @/to, . . ., et la seconde donne 
naissance à une suilo æ, tendant vers l'infini. Ainsi, dans 
la plupart des cas, en écrivant (8) sous la forme (10), on 
n'obtient pas le résultat désiré. 

Nous allons voir quelles sont les propriétés de œ (x) 
qui influent sur la convergence du processus. Comme la con- 
vergence signifie que x, — À — 0 pour x —+ co, il faut 
que |x, — À | diminue lorsque x augmente. Supposons que 
pour x quelconque on ait l'inégalité 


Len X|<Ofrni—X, 8<1. (11) 


Il est alors évident que | x, — X | décroît comme une pro- 
gression géométrique de raison 6, c'est-à-dire qu'il y aconver- 
gence. Substiluons dans le premier membre de (11), 
au lieu de +, et X, respectivement @(x,) et (X). On 
obtient 


[pQu)—p(II<O ana Xf  O<1 (42 
Comme on ne connaît pas la racine X, on ne peut vérifier 
directement la dernière condition et il faut un pou la renfor- 
cer. Ci-dessus nous avons supposé que la racine était localisée 


à l'intérieur d'un certain intervalle. Si pour deux points 
quelconques x’, 2” de cet intervalle la condition suivante 


[p(x)—p()1<O [2-2], 8<1, (13) 


se trouve remplie, la condition (12) se trouve vérifiée à 
fortiori. L'application x = œ (x) satisfaisant à (13) est dite 
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contractante (elle contracte le segment x” — zx’). Il est évi- 
dent que la condition (13) est une condition suffisante pour 
la convergence du processus itératif (9). 

Pour l'estimation de la grandeur © déterminant la 
vitesse de convergence, le plus simple est d'utiliser la 
formule évidente 


8—max|p(x)|, (14) 


où max est pris sur l'intervalle de la localisation de la ra- 
cine. Si l'application x = (x) est contractante, cet inter- 
valle diminuera lorsque #7 augmente, et on peut préciser 
alors la vitesse de convergence. 

Ainsi, lorsque l'on écrit (8) sous la forme (10) tout en 
satisfaisant à la condition (13), on obtient un processus 
itératif convergeant vers la racine. 

Il est évident que la vitesse de convergence détermine la 
qualité de tel ou tel choix de œ (x). En ce sens il va de soi 
que le meilleur est celui pour lequel la grandeur © soit 
minimale. 

Soit æ (x) tel que les équations / (x) = 0 et p (x) — x = 0 
aient une racine commune ÆX. Si cette racine n'est pas 
multiple, /’ (X) 0, au voisinage de cette racine, où la 
fonction / (x) n’a pas d’autres racines, le rapport 

FO E (x) 
f (x) 
est alors une fonction bornée. À chaque fonction (x) corres- 
pond r (x) et inversement chaque r (x) bornée donne naïis- 
sance à 


p (a) = x + r (a) f (x). (15) 


Ce qui nous intéresse, c’est (x) pour laquelle |’ (x) | 
est minimale (en vortu de (14)). En dérivant l'expression 
(15), on a 


q = 1 + r (a) f' (à) + 7 (x) j (a). 
15 


Au voisinage de la racine, la grandeur / (x) est minimalo, 
nous pouvons ainsi négliger le dernier terme et prendre 
l'expression restante égale à zéro. Ceci donne 


ds 
OS FE 
c'est-à-dire, en vertu de (15), la fonction 
\_4+__1(@) 
P (x) = TL f' (x) ’ (16) 


donnant naissance au processus d'itérations, appelé méthode 
de Newton 


En = En — ES. (17) 


La convergence de ce processus est très rapide car 


sp JG) 
p'e)= HE 1 (x) +0 
au fur et à mesure que l’on s'approche do la racine. 

On pout obtenir la formule (17) également d'une autre 
manière (fig. 2), à savoir, à partir d'une certaine approxima- 
tion pour la racine x, 1, lors du calcul de l’approximation 
suivante +, nous allons, au lieu de l'équation f (x) = 0, 
considérer l'équation linéaire 


1 (cn) + fn) (& — tas) = 0 


ne tonant compte que des deux premiers termes du développe- 
ment de la fonction f (x) en sério de Taylor au voisinage du 
point æ,_,. En résolvant cette équation par rapport à x, 
on obtient la formule (17). 

Les avantages de la méthode de Newton apparaissent 
seulement au voisinage do la racine, où œ’ est petit. C'est 
pourquoi dans les calculs il y a lieu, avant tout, de localiser 
la racine par une méthode plus simple et plus grossière, 
puis pour obtenir rapidement une grande précision, utiliser 
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Y=fa-rtfalx-xn. 1) 


Fig. 2 


Ja méthode de Newton. La précision des calculs (nombre de 
décimales retenues) doit être variable, il y a lieu de l'aug- 
menter au fur et à mesure que l'on se rapproche de la racine. 
Comme la transition de +, à +, est simultanément une 
vérification de l'équation æ = ® (x), équivalente à f (x) — 0, 
la précision atteinte peut être estimée d’après le nombre de 
décimales restant inchangées lors de cette transition. 

La méthode itérative de recherche des racines exposée, 
utilisant des applications contractantes, peut sans modifica- 
lions de principe être généralisée à de nombreux problèmes 
plus compliqués. 

Soit à résoudre un système de A7 équations à À7 inconnues 


Jen (20, 20), ,,., a) = O0, m—1, 2, ... A1, (18) 
Désignons par f le vecteur de composantes f(, jf, .., 


..., J80 et par x l'ensemble 2, z%, ,.., æM et écrivons 
le système (18) sous la forme d’une équation vectorielle, soit 


f (x) = 0 (19) 


dont la forme coïncide avec (8). Tout ce qui vient d'être 
dit au sujet de l'équation (8), peut être transposé au systè- 
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me (19), il faut soulement, compte tenu du caractère vecto- 
riel des grandeurs, attribuer un sens nouveau aux désigna- 
tions utilisées. 

Pour l'estimation de la grandeur d'un scalaire on utilise 
son module, pour un vecteur on utilise sa norme. Nous défini- 
rons cetto dernière comme le maximum du module des 
composantes 


Ie] = max | 20. (20) 
m 


Il est évident que si la norme d'un vecteur est égale à zéro 
ceci signifie que toutes ses composantes sont nulles. 


æ = @ (x), (21) 


où æ est le’ vecteur, de composantes qt?, pp, ..., pp. 
Tout comme précédemment la formule (21) donne naissance 
à un processus itératif qui converge si l'application (21) est 


contractante, c'est-à-dire 


lp) p(MI<OIz"—-x 1, O<1. (22) 


Nous allons procéder à l'estimation de la grandeur 6. 
A cet effet considérons la fonction œ‘” (x) sur la droite 
joignant x’ et x”. On obtient une fonction do l'argument 
scalaire s 


D (8) = pin (a' +5 (a"— a°)). (23) 


Les points x’ et +” correspondent aux valeurs s = 0 et 
s = 1. Il est évident que 


PO (a) — pd (x°) = 1p (1) — 1p (0) = ’ (S), (24) 


avec 0ZSs<1. En dérivant (23) par rapport à s, on 
obtient 


, < CLP PENTS CAC 
Ÿ (s)— > oz) (Ce) — (20) 
k—1 


et comple tenu de (24) 
| pr (2°) — pi (21 < 


TS CN max | (eo) — (at) |. 


Compte tenu de Îa PT de la norme (20), on obtient 
M 

Lez") — (a) 11 < Dax MAX > 

Ÿ hf 

Il est cviuent que maintenant c'est la matrice des déri- 


à 1 m) 
gz\4) 


2" — +" |. (25) 


ra ' 
vées — un 
ment qui joue le rôle de œ’. Si l'on définit maintenant la 


norme de cette matrice par la formule 


Ip 1= max D 
mi À 


“du vecteur fonction par rapport au vecteur argu- 


gt") 


a l° (26) 


pour l'estimation de la grandeur © dans l'expression (22), 
en vertu de (25), on peut utiliser l'égalité 


6 = max || (x) |], (27) 
C « 
qui est une généralisation de (14) au cas vectoriel. 


Puis on peut, comme précédemment, construire œ (2) 
à l’aide de (15) 


po =tc+r(x) f(x), (28) 
où r(x) est uno matrice quelconque de fonctions. Pour 
r(x)= — (f(x), (29) 


où (/')-* est la matrice inverse à la matrice des dérivées f” — 
m FL F 

= Er , on obtient la méthode de Newton. Cette méthode 

conserve également ici ses avantages car la matrice qui lui 
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correspond qç’ sera nulle au point de la racine. 11 est facile 
de s'en rendre compte en dérivant (28) et en tenant compte de 
(29). Maintenant, en utilisant la méthode de Newton, au 
lieu de diviser par /” il suffit d'inverser la matrice f’, 
c'est-à-dire lors de chaque itération de résoudre le système 
d'équations linéaires suivant 


(ts) + A (tn) (2 — Tn-s) = 0. 
Lorsque l’ordre du système 17 est élevé, ceci peut être assez 
compliqué. C'est pourquoi généralement on n'utilise cette 
méthode que dans le cas où l’on dispose d'une bonne approxi- 
mation pour les racines, lorsque une ou deux itérations 
donnent une augmentation notable de la précision. 

Il y a lieu de faire la remarque suivante concernant les 
systèmes d'équations linéaires. Si on a besoin de résoudre un 
système du type général il n’y a rien de mieux quelaméthode 
bien connue d'exclusion. Cependant en utilisant les pro- 
priétés spécifiques de certains systèmes particuliers on 
arrive souvent à des méthodes itératives efficaces. 


Problèmes 


1. Construire le processus d'itérations permettant de 
résoudre la racine d'ordre p. 

2, Trouver la rapidité de la convergence de la méthode 
de Newton au voisinage de la racine multiple d'ordre #. 

3. Construire le processus d'itérations de la résolution 
de l'équation ax + b — 0 où 0,1 <a <1, sans utiliser 
l'opération de division. 

4, Soit le système d'équations linéaires Az + b = 0 
de matrice À dont toutes les valeurs propres sont réelles, 
différentes, positives et se trouvent avec certitude à l'inté- 
rieur d’un certain intervalle donné (œ, B). Choisir le nombre r 
de telle sorte que le processus itératif 


Tn = Tnt (Atn-1 + 0) 
converge le plus rapidement possible. 
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$ 7: L'ONCTIONS ET TABLES 


L'une des notions mathématiques fondamentales est celle 
de fonction, celle-ci est utilisée donc dans la formulation de 
la majorité des problèmes. Dans le langage des méthodes 
numériques une fonction ne peut être, de toute évidence, 
représentée que par des suites numériques. Différentes 
méthodes de représentation sont possibles, par exemple, 
à l'aide des cocfficients d'une série par rapport à des fonctions 
ou des valeurs des paramètres dans une formule d'un type 
donné. La représentation la plus courante et la plus univer- 
selle d'une fonction est une table de valeurs. Ainsi, toutes 
les fonctions élémentaires sin x, 1n x, etc., sont des tables. 

Considérons la fonction F (x) et la table correspondante, 
c'est-à-dire deux suites 21, fn (k — 0, 1, 2, ...). Nous 
allons estimer la conformité de la fonction et de la table 
par la précision avec laquelle la table x,, f, permet de 
restituer la valeur de F (x) en un point quelconque x. 

Il est clair que cette précision dépend de la proximité 
de fy à (xx) et de la densité des nœuds de la table des x. 
Le rôle du premier facteur est évident : l'écart de f, à F (x: 
donne la limite de la précision de réproduction qu'il est 
impossible de surpasser sans changer la table. En idéalisant 
le problème nous allons poser 


fn =F (x), k—=0,1,2,... (30) 


Dans ce cas plus la table est détaillée, micux elle décrit 
les particularités du comportement de la fonction F (x), plus 
grande est la précision de restitution, celle-ci dépendant 
évidemment de la méthode do restitution. 

Soit, par exemple, la méthode la plus grossière consistant 
en ce que chaque valeur de f, est attribuée à tout l'inter- 
valle adjacent x, < t << 2x4, (fig. 3), c'est-à-dire que pour 
calculer F (x) on utilise une fonction constante par morceaux 


Po (x) = fn: Th LE That k=0, 1, 2, DE. (31) 
21 


I 
f, | 

f | : 
——. 
( | 

| ] 


Fig. 3 


Pour estimer la valeur de l'erreur P, (x) — F (x) il faut 
avoir en plus de (30) une information supplémentaire sur 
la fonction F (x). Supposons que F (x) soit une fonction lisse 
à dérivée continue. Sur l'intervalle (x,, t,+1) on peut alors 
la représenter comme suit: 


F(x)= EF (er) + FR) @— 2m). 


On voit immédiatement que sur cet intervalle on a 
Po ()—F (x)| max] FI (run), (32) 


c'est-à-dire que l'erreur est de l'ordre de grandeur du pas 
de la table. 

L'interpolation linéaire est une méthode plus précise de 
calcul de F (x) qui consiste en ce qu'au lieu de (31) on cons- 
truit une fonction linéaire par morceaux en utilisant les 
valeurs de gauche f, et de droite /,4, do F sur chaque inter- 
valle (fig. 4). Cette fonction est de la forme 

PAG) = fn (aan), arr (83) 

On peut effectuer l’interpolation par fonctions quadrati- 
ques, cubiques, etc., en utilisant à cot effet trois, quatre, 
etc., points de la table. Considérons le cas général. 
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fut 


Xk-] Xk Xe] X 


Fig. 4 


Soient n + { points 2, 21, ..., +, appelés nœuds de 
l'interpolation et les valeurs correspondantes f,, f1, : . ., fn: 
Formons le polynôme d'ordre nr 


Pa (x) = ù Am (34) 


dont les valeurs aux points x, (i = 0, 1, ..., n) sont 
égales à f;. Posant dans (34) x = x, on obtient un système 
de x + 1 équations pormettant de trouver n + 1 incon- 
nues, à savoir les cocfficients du polynôme 


D am =fr, i=0,1,...,n. (35) 
m=—0 


Lo détorminant do ce système linéaire | 7° | (déterminant de 
Vandermonde) est égal à ll (x; — æj), c'est-à-dire il est 


différent de zéro car tous les x, sont différents. Par consé- 
quent, le système (35) a une solution unique : un ensemble de 
Am- Nous avons démontré qu'il existe un seul polynôme 
d'ordre # (au plus) prenant en (n + 1) points des valeurs 
données. Ce polynôme est appelé polynôme d'interpolation. 
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Il peut s'écrire différemment. La forme (34) est rarement 
utilisée car les expressions des cocfficients a, en fonction 
de æ;, /; sont compliquées. Nous allons écrire le polynôme 
d'interpolation sous la forme de Lagrange, à savoir 


PQ= D ae, (36) 
im 0 


où 
Œ(x)=(r— ro) (tr — 2)... (ris) (ri) (mr — ri). 
Comme q; (xs) = Q pour ik on a de toute évidence 
Ph (tr) = fn. D'un autre coté, chaque g; (x) est un polynôme 
de degré n. Par conséquent leur combinaison linéaire (36) 
cest le polynôme d'interpolation. 

Pour x# — 1 la formule (36) devient 


Pa (2) = fo = + 4 (37) 


To — T1 EE 


c'est la formule d'interpolation linéaire (comparer avec (33)). 

Nous allons estimer la précision avec laquelle le poly- 
nôme d'’interpolation reproduit la fonction F (x). La diffé- 
rence P, (x) — F (x) étant nulle dans les nœuds d'inter- 
polation 2), 21, . - ., tn, le quotient de cette différence 
par la fonction 


q (x) == — à) @ — m1)... (& — 2) (38) 
est une fonclion bornée et l’on peut écrire 
F (x) = P, (x) +R (x) q (à). (39) 


Nous allons estimer Z? (x). À cet cffet considérons la fonc- 
lion auxiliaire 


UV (E) = F(E) — P, (6) — AR (x) q (= CR À) — AR (x)) a (à). 
(40) 


Il est évident que la fonction v (£) s’annule au moins en 
n + 2 points æo, y, : + +, ns 2. Par conséquent, on aura 
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au moins un point & (x) où la (x + 1)-ièmo dérivée de 
v (£) s'annule, si évidemment cette dérivée existe et qu'elle 
soit continue. En dérivant (40) x + 1 fois et en substituant 
— & (x), on obtient: 
0 = FOTO CE (x) —R (x) (n+1)!, 
la (x + 1)-ième dérivée du polynôme P, (£) de degré n étant 
nulle et q (£) étant un polynôme de la forme Et + ,.... 
Ainsi 
5 POP (EC) 
RO ENT 
on obtient ainsi l'expression de l'erreur d'interpolation, 
à savoir 
F (x) P (x) = er FU CE (x)) g (x) (41) 
Le (a+ 1)1 1 

La fonction & (x) reste évidemment indéterminée. 

Si le pas de la table n'est pas supérieur à un certain b, 
c'est-à-dire si 

dau Eh, k=0, 1,...,n—1, 


on a qg(x) — k"*l ol on obtient à partir de (41) 
LFP, (x) | Ken max [FPE (RTE, (42) 
x 


où c, est une certaine constante. On en conclut que pour 
h — 0 l'erreur décroît comme h"*1. 

La relation existant entre la valeur de l'erreur et Île 
degré n du polynôme d'interpolation est plus compliquée. 
Dans la pratique -courante on utilise très rarement des 
formules d'interpolalion de degré assez grand. Ceci pour les 
raisons suivantes. Tout d'abord comme on peut le voir 
à partir de (42), l'augmentation du degré du polynôme 
conduit à une diminution de l'erreur d'interpolation seule- 
ment pour des fonctions très lisses ayant un grand nombre de 
dérivées. Mais généralement nous ne disposons pas d'une 
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information aussi grande sur les propriétés de la fonction 
F (2). Ensuite les valeurs /, sont toujours des valeurs appro- 
chées de F (x,) ne serait-ce que par arrondissement. C'est 
pourquoi les polynômes obtenus à partir de f, et F (x) 
différeront au plus d'une grandeur de l'ordre de f, — F (x). 
De plus les crreurs contenues dans f auront toujours un 
caractère aléatoire, ce qui montre que les fonctions les 
représentant ne sont pas lissos. 

La méthode décrite de restitution de la fonction F (x) 
à l’aide de la table x,, f, à partir du polynôme d'interpola- 
tion n’est pas la seule possible. Nous avons obtenu P, (x) 
(34) à l'aide d’un système do polynômes x” (m — 0, 1, . ..). 
Mais à cotte fin on peut utiliser de nombreux autres systèmes 


m (x). Dans ce cas au licu do (34) il est nécessaire d’'en- 
visager 


D, (x) — 2 mm (t) (43) 


et d'étudier les possibilités et la qualité do l'approxima- 
tion. On peut aller plus loin et trouver la fonction appro- 
ximante sous la forme d'une combinaison non linéaire de 
fonctions de base œp,, (x). Mais il est évident qu'il faut 
avoir pour cela des raisons suffisantes car l'avantage ne pout 
être que dans le rétrécissement du domaine d'application 
possible de la méthode et dans l'utilisation d’une informa- 
tion essenticlle supplémentaire sur F (x) en plus de la table 
des valeurs. 

On peut aborder différemment la question de la restitu- 
tion d'une fonction d’après un ensemble de ses valeurs. Lors 
de la recherche de la fonction d'interpolation nous avons exi- 
gé que ses valeurs dans les nœuds x, coïncident exactement 
avec /,. Mais souvent il suffit d'exiger un écart minimal de 
ces valeurs à celles tabulées. Cette méthode est justifiée par 
exemple dans le cas où / contient à priori des erreurs impor- 
tantes ou si seule la forme do la fonction approximante nous 
intéresse ot non la précision de l'approximation. 
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Nous allons décrire une des méthodes de ce type, à savoir 
la mnéthode des moindres carrés. Soit la table x,, f, (k — 
= 0, 14,...,n). Il y a liou de trouver la fonction 

M 


Par (x) = à AmPm (2), (44) 


où ®@» (x) est un système donné do fonctions (par exemple 
Çm (€) = +”) de telle sorte que la grandeur 


À (Dar (xx) — fa) = 8 (5) 


h=0 


soit minimale. Si M > n, on peut résoudre lo problème en 
trouvant la fonction d'interpolation (44) pour laquelle ô — 
— 0; ce qui nous intéresse, c'est le cas où A7 << n. 

Le choix des coefficients @,, &, . .., a, est arbitraire, 
donc 6 est fonction de ces coefficients. Pour trouver le mini- 
mum de la fonction 6 (@,, &, ..., &sr), annulons les déri- 
vées de celte fonction par rapport à @p, @j, : + «+, Œurs ON 
obtient ainsi le système d'équations 


M n n 
23 @1 23 Pm (tn) Pr (x) = 25 Pm (tn) fs m0, 1, ..., 7. 


En résolvant ce système d'équations linéaires nous trouvons 
les valeurs &;,, &, . .., as déterminant complètement la 
fonction ®,, (x) de (44) qui so trouve être la meilleure 
parmi les fonctions de cette forme pour l'approximation de 
la table x,, f,, si (45) est la mesure de l'écart. 

Nous allons nous arrêter sur certaines applications des 
résultats obtenus. Tello ou telle méthode de correspondance 
entre la table et la fonction permet de réaliser à partir de 
la table différentes opérations fonctionnelles, à savoir l'in- 
tégration et la dérivation. 

Ainsi, s'il faut calculer l'intégrale d'une fonction 
donnée par une table +4, f, (4 — 0,1, ..., Æ), on utilisant 
sur chaque intervalle (x,, x141) la formule d'interpolation 
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linéaire (33) on a 
Th+1 
Pie) de = AE (ans — 2) (46) 
Xh 
En prenant la somme de ces expressions sur tous les inter- 


valles, on obtient la méthode de calcul de l'intégrale. Dans 
le cas où le pas du tableau est constant (x,4, — x, = h)ona 


YK 
| P, (x) dx = (+ fo+ hit ft s…. +frats fx)l (47) 


c'est là la formule quadratique des trapèzes bien connue. 

On peut estimer la précision avec laquelle la formule (47) 
donne la grandeur de l'intégrale de la fonction F (x). Vu (41) 
pour x = À l'erreur d'interpolation est 


e (x) — + F"(E(2))(z— 2x) (t— tn), Ta LT Tu, (48) 


d'où 
Xh+1 
| E (x) dx | œconst max |F”| (trs — xx) *. 
Th EN EX 41 
Th 


En prenant la somme de cette inégalité sur tous les interval- 
les compte tenu de Xh — xx — x,, on obtient 


X 
|| e(r) dæ|<const max |7"|h?, (49) 
. XOSXEX je 


c'est-à-dire la précision de la formule des trapèzes (47) 
est de l'ordre de h?. En utilisant d’autres fonctions d'inter- 
polation P, (x), P, (x), etc., on obtient les formules quadra- 
tiques des rectangles, de Simpson, etc. 
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Le calcul de la dérivée d'une fonclion tabulée est égale- 
ment très simple. En utilisant par exemple l'interpolation 
linéaire (33), on obtient 

AE Ps frs fn 
de dx That Th 
En dérivant l'expression de l'erreur d’interpolation (48) on a 


; ThLTETh+1e (50) 


De 7 A7 (€) (r)) E' (a) (ee) (en) 
+ + FE" (E(x)) (27 — 23 — Zn). 


Le second terme dans le second membre est de l'ordre de 
Lt — Th C'est-à-dire de k. C'est la précision du calcul 
de la dérivée assurée par la formule (50), à l'exclusion du 
point central de l'intervalle & = (x, + 2,+,)/2. En ce point 
le second terme s’annule et par conséquent la formule (50) 
donne la valeur do la dérivéo en ce point à k? près. 

Lorsque l’on essaie de calculer la dérivée seconde par 
interpolation linéaire, on obtient: 


—— PO) mt 

dr? dx? 
ce qui ne convient évidemment pas. Ceci se trouve en accord 
avec le fait que 


= FE (æ))+ 


c’est-à-dire que l'erreur est finie, ne décroît pas lorsque 
diminue. Pour le calcul des dérivées d'ordre plus élevé il y 
a lieu d'utiliser une interpolation d'ordre plus élevé. 

Nous avons envisagé seulement le cas des fonctions d'une 
seule variable. Lorsque l’on passe à des problèmes multi- 
dimensionnels, les principes de base des méthodes exposées 
restent valables, mais de nombreux autres problèmes appa- 
raissent. 
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Problèmes 


1. Trouver c, en fonction de 7 dans la formule (42). 
2. Soit la table 


Trouver l’approximation de cette table par une fonction liné- 
aire en utilisant la méthode des moindres carrés. Puis trouver 
son approximation par une fonction linéaire en utilisant 
pour l'estimation de l'écart non pas (45) mais 


max | D (x») — fa | = 6. 
h 


Comparer les fonctions obtenues avec le polynôme d'inter- 
polation du second degré obtenu à partir de ces mêmes 
points. 

3. Trouver les formules quadratiques en utilisant les 
polynômes d'interpolation de degrés zéro et deux. Trouver 
l'estimation aussi précise que possible de l’ordre de gran- 
deur de l'erreur de ces formules. 

4. Trouver la formule de calcul de la dérivée seconde 
d'une fonction tabulée. Estimer l'erreur. 

5. Une fonction de deux variables F (x, y) est donnée 
par la table suivante 


fo, 1 | J-1,0 | fo,-1 


Jrouver la formule quadralique pour le calcul de l'inté- 
grale 


| | F(x, y) dx dy, 
x2+y2<h9 


en utilisant dans chaque octante l'interpolation linéaire sur 
deux variables. 


3 ÉQUATIONS DIFFÉRENTIELLES 
$ °  ORDINAIRES 


Nous allons étudier un problème typique, bien que pas 
très compliqué, pour les équations de ce genre. Il y a lieu 
de trouver la fonction U (f) satisfaisant pour t => 0 à l'équa- 
tion 


au , 
SAR "(E, U) (91) 
prenant pour £ = 0 une valeur donnée 
U (0) = U, (52) 


En théorie des équations différentielles ordinaires, si 
le second membre de (51), F'(t, U) satisfait, en tant que 
fonction de ses arguments, à certaines conditions de régu- 
larité, la solution du problème (51), (52), U (t) existe, elle 
est unique et se trouve être une fonction régulière. Nous 
allons supposer que ces conditions soient vérifiées. 

Les cas où U ({) peut s'exprimer à l'aide de fonctions 
élémentaires ou par des intégrales de ces fonctions sont 
exceptionnels. En général, seules les méthodes numériques 
permettent de résoudre le problème (51), (52). Elles donnent 
évidemment une information limitée, approchée sur la solu- 
tion mais sont en revanche universelles. Nous allons exposer 


ici la plus simple de ces méthodes, à savoir la rmnéthode 
d'Euler. 


31 


On remplace avant tout le domaine de la variation conti- 
nue de l'argument { > 0 par un ensemble discret de points 


{= NT, 0:12 ue (53) 


où t est un certain nombre donné petit, dit le paramètre de la 
méthode numérique. Au lieu de la fonction U (t) nous allons 
considérer la table 


bot MU 2e. | (54) 


Comme par définition de la dérivée dU/dt est la limite du 
rapport (U (t + t) — U (t))/t pour t —+ 0, en remplaçant 
la dérivée par ce rapport fini, on obtient au lieu de l'équa- 
tion différentielle (51) l'équation aux différences 

nn EF, un), n=0,1,2,... (55) 


T 
ou 

Unti = Un TE (tn, Un); n—0, 1, 2,..., (55) 

Ainsi, compte tenu de (52) et posant 
Ugo = U; (56) 
on peut à partir de (55”) trouver successivement tous les uw. 

On a ainsi trouvé l'algorithme de calcul. 
Naturellement, maintenant l'algorithme est d'autant 
meilleur que le tableau f,, #, reproduit mieux la solution 


exacte du problème initial (51), (52), à savoir la fonction 
U (t). Pour le voir, posons 


ln - U ({n) “a Oltns (57) 


et essayons d'estimer la grandeur de l'erreur ôuw,. 
En substituant (57) dans l'équation aux différences (55), 
on obtient 


dUnst —OUn = F(t,, Ù (En) + Bt.) — En On), (58) 


T 
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Nous allons estimer le deuxième membre. Comme U (t) 
est une fonction régulière et qu'elle satisfait à (51), on a 


U (tn41) =U (x) +7 (F). +0, (7?) = 
=U(t)+7TE (tft, U (t:)) +0 (x), 
d'où 


Una) 7 On) = F(t,; U(tr)) +0 (x). (59) 


En comparant cette égalité avec (55), notons que la solu- 
tion exacte de l'équation différentielle du problèmeinitial 
satisfait à l'équation aux différences (55) à © (rt) près. 
En substituant maintenant (59) dans le deuxième membre 
de (58), on peut écrire 


Mn Un (fn, U (tn) + Gun) —F (tn, Ü (tn}) + O (r). (60) 
Vu les condilions de régularité imposées à Æ (t, U), on a 
l'inégalité 

LFGn, Ut) +8u,) — Ft, U())IS M] 6u,], (61) 


où À est une constante. En utilisant (61) on obtient à partir 
de (60) 


| Étn 41 [SL + TA2) | ôu, | + O (r?). (62) 


La formule (62) donne l'accroissement de l'erreur pour 
un pas. Donc, pour l'erreur accumulée durant N pas, on a 


ôux|<O (T°) + (+ TA) | ôux|< 
LOT) + (1 + T42) (0 (x) + (1 + TA) | Eux 2 & 


(+ TANT O (6) + (1 + T1)" | wo |, 
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d'où, en pronant la somme de la progression géométrique 
à raison À + TA, on a 


14 Tu)" —1 
un 0 2) + (4 + TM) Ouol. (63) 


Ce qui nous intéresse c'est la valeur de l'erreur pour 
des v petits pour { donné quelconque. Posant £{ = {Vt et 
tenan. compte de ce que 


A+iM)" = (A + rat)" et pour te 0, 
on obtient à partir de (63) 


| du (#)|SeMtO (x) + et] Eu |. (64) 
Comme ôus = uo — Uy = 0, il en résulte que pour + —+ 0 
l'erreur ôu sur un intervalle fini quelconque tend vers zéro, 

Nous avons ainsi démontré que pour t suffisamment petit 
la table obtenue par la méthode d’Euler donne une approxi- 
mation aussi précise que l'on veut do la solution du problème 
initial (51), (52). 

Nous n'avons pas obtenu le critère permettant de dire 
quel est + que l’on pout considérer « suffisamment petit » 
pour obtenir une précision donnée. Nous n'avons que simple- 
ment établi la convergence de la solution approchée vers la 
solution exacte pour Tt —+ 0. 

On aurait pu en utilisant au lieu des désignations O (x), 
O (tr?) des expressions plus détaillées donner une estimation 
quantitative de la grandeur de l'erreur. A cet effet on 
n'aurait eu qu’à surmonter des dificultés purement techni- 
ques. Mais une telle estimation serait très surélevée et ne 
saurait être utilisée pour calculer l'erreur exacte. 

C'est pourquoi, en cas de besoin, pour trouver la préci- 
sion accessible, on procède comme suit. Comme & — U pour 
tT—+ 0, à partir d’un certain + les premières décimales signi- 
ficatives u cessent de changer, elles correspondent donc à la 
solution exacte. En effectuant les calculs avec différents t 
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et en comparant les résultats obtenus, on peul so faire une 
idée de la précision obtenue. 

De ce point de vue les méthodes accusant une forte dépen- 
dance de la valeur de l'erreur du paramètre t, i.e. les métho- 
des de haute précision, sont plus commodes. La méthode 
d'Euler a une précision d'ordre minimal, u — U = O0 (x). 
Il est évident que ceci provient de l’approximation grossière 
employée consistant à remplacer l'équation différentielle 
par une équation aux différonces. 

Pour estimer la qualité de l'approximation on détermine 
la précision avec laquelle la solution du problème initial 
satisfait à l'équation aux différences. En comparant (55) et 
(59), on voit que dans le cas présent l’approximation est de 
l'ordre de © (x). Cette méthode d'estimation de l’ordre de 
l'approximation est quelque peu arbitraire. Par exemple si 
au lieu de (55) on utilise (55°), l'ordre de l’approximation 
sera © (ti), vu la norme différente de l'équation aux différen- 
ces (multiplication par t). Nous allons utiliser la norme 
naturelle de l’équation aux différences pour laquelle cette 
dernière devient, à la limite, l’équation différentielle ini- 
tiale. Ainsi, (55) (plus exactement (59)), pour t —+ 0, devient 
(51), alors que (55°) devient l'égalité U = U où les pro- 
priétés du problème initial n'apparaissent pas du tout. 

Considérons l'équation aux différences 


n+i— ‘°n 1 
RAR EF (bn, Un) HF (tntis Unts)) (65) 
et eslimons l'ordre de l’approximation de l'équation diffé- 
rentielle (51) correspondante. En substituant dans (65) au 
lieu de w,, uh+, les valeurs U (£,), U (t1+1) et en utilisant 
les développements 
du 2 fd2U 
UÙ(tn41) = U (tn) +7 (T).++ Gr), +06), 


PF (bnras U (En) = F (bn D (bn) +7 (Sr), +0 (2) 
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on obtient l'égalité 


(+ T Te) x dF 


— (F+-5 7), +00). 


UT ni 


dt 2 dt? 


Comme U (t) satisfait à l'équation différentielle (51), ainsi 
qu'à 


d2U __dF 
dtè — dt 


en découlant, on en conclut que pour & = U (65) est satisfait 
à O (1°) près. 

Ainsi l'équation aux différences (65) est une approxima- 
tion de l'équation initiale (51) à © (x°) près. Le fait que 
la solution conserve cette même précision, c'est-à-dire que 
u — U = 0 (1°) n'est pas évident, et il y a lieu de démontrer 
celte propriété. Mais nous n'allons pas nous arrêter sur ce 
point. 

À la différence de (55), la formule aux différences (65) 
ne permel pas en général d'exprimer d’une manière explicite 
Un+1 en fonction de u,. C'est une équation par rapport à u,,+1. 
Pour trouver la solution, on peut utiliser tel ou tel proces- 
sus d’itérations, d'autant plus que l'on a toujours une bonne 
approximation initiale u,. Cependant il est inutile d'essayer 
de trouver la valeur exacte de #,+1, ce serait pure perte de 
temps, car l'équation elle-même est donnée à © (1?) près. 
C'est pourquoi on peut se limiter aux deux itérations sui- 
vantes. Tout d'abord on calcule la première approximation 


Ut par la formule de la méthode d'Euler 
Unki= Un + TF (tn, Un) (66) 


que l’on substitue ensuite dans le deuxième membre de (65). 
On trouve ainsi la valeur définitive de u,+1, 


Uni = Un + CE (Enr Un) + F (ntir Unis). (67) 
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En fait, ceci signifie qu'au lieu de (65) on utilise l'équa- 
tion aux différences 


East EF (En, Un) + F(tnr1s Un + TF (ln, Un))). (68) 


T 


ll est facile de voir que celte équation, tout comme (65), 
est une approximation de l'équation initiale (51) à © (x°) 
près, donnant l'expression explicite de w,+, en fonction 
de «,. 

Remarquons que lorsque l'on effectue le calcul à l'aide 
de (66), (67), on peut contrôler la précision du résultat 


obtenu en comparant les valeurs #,41 et u,41 sur chaque pas 
de calcul. Une trop grande ou une trop petite valeur de la 
différence un+1 — Una Signifie qu'il y a licu de diminuer 
ou d'augmenter le pas t. En changeant le pas dans le sens 
convenable, on peut maintenir la précision à un niveau 
donné. 

En généralisant la méthode exposée de l'équation aux 
différences on peut obtenir des méthodes dont la précision 
est d'un ordre encore plus élevé. Si l'on prend en considéra- 
tion plusieurs valeurs intermédiaires de « calculées succes- 
sivement, on arrive à des méthodes du type de celle de Runge- 
Kutta (voir problème 2). Une autre généralisation consiste 
à utiliser pour obtenir w,+, non seulement w, mais également 
les valeurs connues w, 1, Un», . . . Le principe général de 
ces méthodes (du type d’Adams) est le suivant. À partir de 
la suite donnée F (t,, un), F (fi-5s Un), + «+ F'(tn-ns Un -x) 
on trouve le polynôme d’'interpolation P, (t). En l'utilisant 
au lieu de F sur l'intervalle #,, {+1 (extrapolation), on peut 
écrire en intégrant (51) 


fn+1 
Un+ti — Un = Ph (£) de. 


în 
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Comme P, (4) cest une fonction linéaire des valours do F aux 
points mentionnés, après intégration on obtient 


R 
Unti= Un +T 2 CF (tn-t, Un-1); (69) 


c'est-à-dire une formule aux différences. Pour k — O0 (69) 
devient la formule de la méthode d'Euler (55). Connaissant 
Un+1 On peut, comme dans la méthode d'Euler, introduire 
une correction en ulilisant le polynôme d'interpolation 
tenant compte de la valeur w,4, qui vient d'être obtenue. 
Cetto méthode permet do réaliser des méthodes numériques 
d'une précision quelconque et ceci d'une manière très écono- 
mique car l'augmentation de la précision s'obtient en utili- 
sant les valeurs de z dont on dispose déjà. 

Néanmoins ces méthodes, très répandues lorsque les 
calculs étaient effectués à la main, ne sont presque pas 
utilisées à l'heure actuelle. Ceci pour des raisons très carac- 
téristiques pour les méthodes modernes de calcul numérique. 
En effet, on peut utiliser la formule (69) seulement à partir 
d'un certain # — k. Par conséquent les valours 41, u2, . .. 

.., tn doivent être obtenues différemment, par des for- 
mules spéciales. Des difficultés analogues apparaissent 
lorsque l'on change le pas de l'intégration. La nécessité 
d'étendre l'algorithme pour tenir compte de plusieurs cas 
exceptionnels rend ces méthodes peu utilisées. 

Tous ce qui vient d'être dit ci-dessus peut, tout natu- 
rellement, être généralisé à un système d'équations différen- 
tielles (et par conséquent à des équations d'ordre plus élevé). 
En particulier, si U, F, u ne sont pas des grandeurs scalaires 
mais des grandeurs vectorielles (51) devient un système 
d'équations différentielles, et les formules aux différences 
(55) et autres décrivent les méthodes d'intégration de ce 
systèmo. Les résultats de l'étude de la convergence de l'ap- 
proximation restent également bien que l'étude elle-même 
se trouve être un pou plus compliquée. 
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Lors de l'élaboration des méthodes numériques nous 
avons systématiquement et essentiellement supposé la fonc- 
tion F (t, U) régulière. Si l’on renonce à cette hypothèse et 
que l’on suppose que la fonction F ({, U) présente des singu- 
larités, on ne sera pas de toute évidence obligé d'élaborer une 
méthode spéciale de solution qu'au voisinage de ces singula- 
rités. Par exemple, si la fonction F (t, U) présente un point 
singulier, du type indétermination 0/0, il y a lieu d'isoler 
les termes principaux ? ({, U), de trouver la solution analy- 
tique approchée au voisinage do ce point et de l'utiliser pour 
passer le point singulier, on juxtaposant à la table #,, u, 
obtenue par la méthode numérique. 

Enfin, nous allons nous arrêter sur le rôle des erreurs 
d'arrondissement. Comme nous l'avons déjà noté, toute 
équation aux différences porte une erreur d'approximation. 
Si l'imprécision apparaissant sur un pas par suite de l'arron- 
dissement ne surpasse pas la précision de l'approximation, 
cette relation reste vraie pour les erreurs accumulées dues 
à l’une et l'autre de ces causes. Par conséquent, si avec 
la diminution de % on augmente la précision des calculs, 
la convergence se conserve. Considérons par exemple la 
méthode d'Euler. En fait, le processus de calcul nous donne 
non pas une suite w, satisfaisant à (55°) mais une certaine 
suite voisine w,. Cette dernière est obtenue par la formule 
suivante 


Un+i — Un +TE (Es, Un) + Ôns 
où 6, est l'erreur admise lors du calcul du second membre. 
Elle est due à l’imprécision du calcul de F (£,, u.) et à 
l'imprécision do la multiplication par +. Il est évident que 
6, caractérise la précision des calculs. Si 6, = 0 (T°), alors 


4 


u, satisfait à la même équation que U (4), c'est-à-dire 


à (59). Par conséquent, #, — u, = O (xt) ot le processus 
de calcul réel mené avec une précision excédentaire Ô — 
—= 0 (x?) pormet d'obtenir la solution à © (x) près. 
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Problèmes 

1. Montrer que la solution obtenue par la méthode 
d'Euler avec les corrections ((66), (67)) coïncide avec la 
solution exacte à O (x?) près. 

2. La méthode de Runge-Kutta d'intégration de l'équa- 
tion (51) est décrite par les formules suivantes: 


Un —= Un ++ (/: + 2/2 + 2/3 + fe), 
où 
J ES F (Es Un); 


f2=F (a+, ui+T fs) , 
Ja=F (++, un +z fa) ' 
li =F(t4 +T, Un + Tfa)e 
Trouver l'ordre de l’approximation de cette méthode. 
3. Montrer que la méthode d'Euler pour un .système 
linéaire 
I 
Fire ayjU;, ES À 
j=i 
donne la solution à © (rt) près. 


4, Décrire l'algorithme de calcul donnant la solution 
du systèmo 


dx x 

dy 

do re 
satisfaisant aux données initiales { — 0, x — 0, y — 0 et 
prenant pour un certain { — 7' des valeurs données x = X, 


y = YŸ. La fonction f ({) est régulière et donnée par une 
table admettant une interpolation linéaire. 
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5. Pour la solution du problème 
SP +MU=0, U(0)=1 
utiliser les deux équations aux différences 
HR + Mun = 0 
et 


tu — 
ne + Mu, =0. 


Estimer l'ordre de l'approximation dans les deux cas. 
Comparer les solutions approchées, obtenues pour 7 fini, 
avec celle exacte. Lequel des deux cas est préférable pour 
M grand lorsque seul le caractère qualitatif de la solution 


exacte nous intéresse ? 


Chapitre Il 


ÉQUATIONS AUX DÉRIVÉES 
$ À, PARTIELLES 


Tout ce qui va suivro sera consacré aux méthodes numé- 
riques de solution des équations différentielles aux déri- 
vées partielles. Lorsque l'on passe d'une à deux 
et plus variables indépendantes, la diversité et la complexité 
des problèmes augmentent brusquement. N'ayant pas la 
possibilité de voir on détail tous les aspects et les méthodes 
de cette théorie en plein essor nous nous limiterons à l'étude 
de certaines questions de principe. 

Commençons par le problème le plus simple. Dans le 
domaine 


—o0 <Lr<o, t1>0 (70) 
il y a lieu de trouver la fonction U (x, t) satisfaisant pour 
t> 0 à l'équation différentielle 
+ a = F (x, t) (71) 
et prenant pour & = 0 des valeurs données 
U (x, 0) = ® (x). (72) 


La solution exactè de ce problème est donnée par la 
formule 
t 
U(x, t)=O(z—at)+ | F(æ—at+at,t}dt, (73) 
0 
comme on peut facilement s'en rendre compte. Mais pour le 
moment ce problème ne nous intéresse qu'à titre d'exemple 
élémentaire, permettant de montrer la plupart des propriétés 
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des méthodes numériques d'intégration des équations dif- 
férentielles aux dérivées partielles. Ci-dessous nous y re- 
viendrons encore une fois. 

Pour l'élaboration d'une méthode numérique de résolu- 
tion du problème (70), (71), (72) il faut avant tout remplacer 
le domaine des varialions continues des arguments (70) par 
un réseau de calcul (fig. 5), c'est-à-dire par un ensemble 
discre des points de coordonnées 


Zn = kk, k — 0, +1, +2, 0) 
NT, = 0: 1:23 (74) 
On voit ainsi apparaître deux paramètres v et k, ce sont les 
pas du réseau. Au lieu des fonctions U (x, t}, F (x, t), D (x) 
on va envisager les fonctions discrètes, c'est-à-dire des suites 
numériques u%, fh, ®x Correspondant aux points du réseau x,, 
{" (74). Puis, en remplaçant los dérivées partielles figurant 
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dans (71) par les rapports aux différences on a 


n+i n n n 
lu — ll U — 
k k h+i k n 
= A — = fh (75) 


pour chaque paire d'indices k, n, c'est-à-dire pour chaque 
point du réseau de calcul (74). Enfin, au lieu de (72) on écrit 


ui = Pr: (76) 


Ainsi nous avons remplacé le problème (70), (71), (72) 
par le problème numérique aux différences (74), (75), (76). 
Il est évident que la méthode utilisée n'est pas unique. Le 
nombre de possibilités se présentant pour la résolution des 
équations aux dérivées partielles est bien plus grand que dans 
le cas des équations différentielles ordinaires. Mais pour le 
moment nous en resterons là. 

L’algorithme de calcul donnant les valeurs de uf est 


très simple. Nous allons résoudre (75) par rapport à uf%*1, 


ut} = (1 + a +) Uh — + uysi + Tfr (77) 


Si les grandeurs u* pour une certaine couche »# de points sont 
données la formule (77) permet de les calculer pour la 
(x + 1)-ième couche de points. Comme uf sont donnés, ils 
permettent de calculer «uk, puis uf et ainsi de suite. 

Nous allons passer maintenant à la question essentielle, 
à savoir à quel point la solution uw obtenue par une méthode 
numérique est voisine de la solution exacte du problème 
initial U (x, t). Il est évident que ceci ne peut avoir lieu 
que si t et À sont petits. 

Posons 


uh = U (xp, €") + Oui (78) 
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et substituons cette expression dans la formule aux diffé- 
rences (75). On obtient 


un +1 Eu . ôuY 4 — Ôuy 


6 d h 
Dutt—UR Uri Un 
= (+ HE), 9) 


Ici et plus loin UX£ désigne U (x,, t"). Nous allons estimer lo 
deuxième membre de (79). En supposant que U (x, t) est une 
fonction régulière, pour + et À petits on peut écrire 


URt' = Ur + (5) +0), 


(80) 
Ua = Un + h (5) + 0 (ue), 
d'où 
n+i _ym n 
ON A ae (+ La a) +O(r, h). 


Comme UÙ (x,t) satisfait à (71) ot F'(x,, 1") = fx, la der- 
nière égalité signifie que 
Un+I_ Ur U} H= 


” La — fr +O (x, h). (81) 


En comparant (81) avec . on voit que la solution du 
problème initial U (+, ?) satisfait à l'équation aux différen- 
ces (75) à O (t,h) près, c'est-à-dire que nous avons une 
approximation. 

En substituant (81) dans le deuxième membre de (79), 
on obtient l'équation déterminant ôuw 


Gun + — ôuy (THE — ou} 
mr 2 À — 
T h 


—O(x, h) (82) 
Ou 


Su +! — (1+a+ ) Suk — à + Euh 41 + TO (x, h). (83) 
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Envisagcons tout d'abord le cas où a, t, À satisfont aux 
inégalités 
OZ — a + <1. (84) 


Dans ce cas les cocfficicents près de ôu% ct ôux4, dans le 
deuxième membre de (83) sont positifs, et l’on peut écrire 


| Guy +! |S (1 La +) | ôux | + — +) ôuxsi [+ TO (x, h)< 
<max(|ôux|, |ôux41l) + TO (Tr, h). 
Introduisons la désignation 


| Su" || = max | Our |, (85) 
k 


l'inégalité précédente donne alors 
[ ôu+t [<< I] Bu” + TO (x, Xe), (86) 


c'est-à-dire que l'écart maximal ô4 pour un pas augmente 
de TO (t, À) au plus. Pour W pas on aura 


[ôu" 1II êu° 11+ NrO (x, k). (87) 


Fixons { — Nr fini quelconque et faisons tendre +, À vers 
zéro, N tendra alors à l'infini. Comme ôuÿ — D, — y, on 
obtient à partir de (87) 


[| ôu (t) |] = O (t,.h). (88) 


Nous avons ainsi démontré que si, lorsque +, À tendent 
vers zéro, les conditions (84) se trouvent vérifiées, la solu- 
tion du problème aux différences (74), (75), (76) tend vers 
la solution du problème initial (70), (71), (72). 

Considérons maintenant lo cas opposé lorsque +, À ten- 
dent vers zéro de telle sorte qu'au moins l'une des condi- 
tions (84) ne soit pas vérifiée. Il se trouve que dans ce cas, 
en général, il n'y a pas de convergence. Nous allons montrer 
ceci par le raisonnement simple suivant. 
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Fig. 6 


Déterminons pour u” le domaine de dépendance. Puisque 
uN s'exprime en fonction de fN"1, uN71, uN-1 qui s'expriment 
à leur tour en fonction de f}"3, fN-3, uN-3, uN-3, N°3, etc. 
(fig. 6), en prolongeant co processus on peut exclure toutes 
les valeurs intermédiaires ux et oxprimer uŸ directement 
en fonction do f% et p,. Les points X, n, correspondant aux 
valeurs f ot o utilisées forment de toute évidence un triangle 
représenté fig. 6, qui se trouve être le domaine de dépen- 
dance de uŸ. Ses côtés sont formés des segments de droites 


re _ bp? _ 
z=0, t=0, tee, (89) 


c'est-à-dire qu'il est determiné par deux paramètres x 
et t/h. Si ces paramètres sont constants pour Tt, À —+ 0, alors 
le domaine de dépendance reste inchangé. 
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Nous allons maintenant déterminer le domaine de dépen- 
dance de la solution exacte en ce même point, soit U (0, Mr). 
Comme on peut le voir à partir de (73), U (0, Vx) se trouvo 
entièrement déterminé par les valeurs F (x, t) sur la droite 


xz=a{(t — Nr), 0O<t< NT, (90) 


et par la valeur ® au point d'intersection de cette droite 
avec l'axe &. 

Supposons que la droite (90) se trouve à l'extérieur du 
lriangle (89). Dans ce cas U (0, Nr) et uN seront détermi- 
nés par différents facteurs indépendants et nous n'avons 
aucune raison pour penser qu'il seront voisins. ln changeant 
par exemple les valeurs / el ® seulement au voisinage de 
la droite (90) nous obliendrons différents U (0, Nr). 
De plus uŸ ne sera pas influencé par ces variations. 

Trouvons les conditions d'appartenance de la droite (90) 
au triangle (89). Pour { << Tr donné, on a pour les points 
intérieurs du triangle (89) 


t 
Ogrsh(N——). 
En substituant au lieu de x son expression (90) on obtient 
Oa(t—NT)&E(NT—0) 


ou en simplifiant par NT — 1! 


0< as, 


ce qui de toute évidence est identique à (84). 

Ainsi nous avons montré que lorsque les conditions (84) 
ne sont pas vérifiées, en général, il n’y a pas de convergence. 

[1 y a lieu de faire la remarque importante suivante. 
Le problème aux différences (74), (75), (76) donne une ap- 
proximation de l'équation différentielle initiale (70), (71), 
(72) que la condition (84) soit vérifiée ou non. Cependant il 
se trouve qu’une seule approximation ne suffit pas pour la 
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convergence de u vers la solution exacte U (x, ?). La condi- 
tion (84) est, dans le problème présent, la condition supplé- 
mentaire assurant la convergence. Lorsque cette condition 
n'est pas vérifiée, l'écart de ux à U (x, t) peut être quel- 
conque. 

Néanmoins il est intéressant de voir le caractère de la 
solulion u} obtenue pour +, À donnés et finis et ce qui se 
passe lorsque +, À diminuent. Pour s’en rendre compte sans 
effectuer de calculs nous allons procéder comme suit. Reve- 
nons à la formule (83) décrivant le processus d'évolution 
de l'erreur d'une couche à l’autre. Vu la présence de tO (x, b) 
dans le deuxième membre, l'ordre de l'erreur n'est pas 
inférieur à cette grandeur. Quant à l'erreur ôu#, c'est une 
certaine fonction compliquée de l'indice 4. Supposons qu'on 
puisse l'écrire sous la forme d une somme dont l'un des 
termes est de la forme e (—1)* où de toute évidence e est 
de l'ordre de TO (x, k). Nous allons voir comment change 
cette seule composante de l'erreur, c'est-à-dire posons 


ôuy —e(—1)* (91) 


et calculons ôuft", ôux'?, . .. En idéalisant quelque peu le 
problème, nous ne tiendrons pas compte de l'influence du 
terme TO (rt, h) sur les pas suivants. On obtient 


Gui = (1+a)e(—1}#—a +e(—1} = 
= (1420 +)e(—1y, 


autrement dit, la fonction de la forme (91) sur la couche sui- 
vante conserve sa valeur multipliée par 1 +- 2at/h. Il est 
évident qu'après # couches ce facteur est égal à (1 + 2at/h)\, 
soit 


sui" = (1424) e(—1ÿ, (92) 
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et par conséquent l'évolution de la composante étudiée de 
l'erreur est donnée par la grandeur 1 + 2at/h. Si 


1424 + |<t, (93) 


l'erreur décroît, dans le cas contraire elle croît exponen- 
tiellement. 

Ainsi la solution uf contenant cette erreur perd rapide- 
ment son sens devenant une suite chaotique de termes très 
grands. Même si e est petit, cela ne change rien et ne peut 
que reculer l'instant de la catastrophe. L'effet décrit est 
appelé instabilité. 

Il est facile de voir que (93) est de nouveau la même 
condition (84). Mais ceci signifie que l'absence de conver- 
gence et l'instabilité sont dues à une même cause. En utili- 
sant (92) pour un segment fini { — + on obtient pour le 
cas où la condition (93), ou (84), n'est pas vérifiée: 


|ôux+! =] 1 + 2a +["+0 (t, h)— oo pour T, k +0. 


Le calcul suivant le schéma aux différences instable non 
seulement ne donne pas un résultat voisin de la solution 
exacte mais est impossible. Dans ce cas la diminution de t,A 
no fait qu'accentuer cette instabilité. 


Problèmes 


1. Sia > 0, les conditions (84) ne sont jamais vérifiées 
quels que soient t, k. Comment faut-il changer la formule 
aux différences (75) dans ce cas? 

2. Elaborer et étudier la méthode aux différences de 
solution du problème (70), (71), (72) dans le cas général de la 
variable a, «a — a (x, t). 

3. Pour résoudre le problème 

2 
Ter Ur 0)=D(x) 


ôt Or? 
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envisager l'équation aux différonces 


n+1 n n n n 
nt —Uy Uhy—2u; +uy_; 0 
= QR 


Etudier son approximation (pour 2t/h? & 1) et sa stabilité 
(sur la fonction e (—1)"). 


$ D: APPROXIMATION ET STABILITÉ 


Le schéma de principe de l'étude de la convergence 
faite ci-dessus est typique pour un grand nombre de pro- 
blèmes et peut être décrit de la manière suivante. Supposons 
qu'un certain problème de base différentiel soit remplacé 
par un problème aux différences du type 


lu = f. (94) 


u, / sont ici des fonctions discrètes u%, fà et l'un opérateur 


linéaire aux différences dépendant des paramètres +, À. 
En particulier le problème envisagé (75), (76) s'écrit sous la 
forme (94) à condition de poser 


un+i —u} 44 —Uf 
lu = T ho? (95) 
0 
Uk 
et 
k 
L 96 
EE (96) 


Pour l'étude de la convergence de la solution du problème 
(94) à la solution du problème initial, à savoir à la fonction U, 
on pose 


u = U + Ôôu 
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et en substituant cette expression dans (94) on obtient 
l'équation détorminant ôu 


lôu = f — LU. (97) 


Il est évident que la convergence, c'est-à-dire la tendan- 
ce de Ôu vers O0, aura lieu si, premièrement, par le choix des 
paramètres +T, À le second membre de (97) peut être rendu 
aussi petit que l'on veut, et deuxièmement, s'il reste petit 
pour la solution de l'équation (97), c'est-à-dire pour la 
fonction discrète ôu. 

La première de ces conditions 


IU—f—+0 pour t,h—+0 (98) 


est la condition d'approxzimation que nous connaissons déjà 
(comparer avec (81)). Elle caractérise la relation existant 
entre le problème aux différences et le problème différen- 
tiel et établit le voisinage de ces problèmes. 

La réalisation de la seconde condition dépend seulement 
des propriétés du problème aux différences, à savoir : l’opéra- 
teur aux différences / doit être tel que pour +, quelconques 
la solution du problème (97) soit du même ordre de gran- 
deur que lo deuxième membre, c'est-à-dire 


ôu — f — IU. (99) 


Peu d'opérateurs aux différences sont doués de ces pro- 
priétés. Comme nous l'avons montré ci-dessus pour l'opé- 
raleur / (95), la condition (99) est vraie si a, t, k satisfont 
à la relation (84). Dans le cas contraire la condition (99) 
n'est pas vérifiée. Comme dans ce cas il y a instabilité, en 
attribuant à ce terme un sens plus large, la condition (99) 
est appelée condition de stabilité. 

Remarquons que les équations (97) et (94) ne diffèrent 
que par les désignations: Ôu, f — IU au lieu de «, f. La con- 
dition de stabilité (99) peut donc s'écrire sous la forme 


u — f, (100) 
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et la solution du problème (94) doit satisfaire à cette condi- 
tion. Soulignons que l'opérateur ! dépend des paramètres 
T, hot la stabilité signifie que la relation (100) cest vérifiéo 
pour des t, À aussi petits que l'on veut. 

Pour donner un sens exact aux relations (98), (99), (100) 
il y a lieu de mentionner le mode d'estimation des grandeurs 
y figurant «, f, {U, elc., c'est-à-dire d'introduire une norme 
de ces fonctions de réseau [lu ||, [[/ 1], I ZU I, etc. Dans 
l'exemple envisagé ci-dessus nous avons utilisé en qualité 
de norme le maximum du module des valeurs de la fonction, 
autrement dit, l'égalité 


ôu = O (t,h) 
signifie quo 


[ôu|| —max|ôurx|—0O(t, k). 
h,n 


D'autres définitions de la norme sont également possibles, ce 
qui importe, c'est qu'ello nous satisfasse en tant que méthode 
de mesure de la valeur exacte de la fonction de réseau. 
Ainsi la norme || ôu || — max k | ôux | ne convient pas car, 


1 nn 
dans ce cas, || ôu [| — O0 pour À —+ O0 même si ôu; restent 
finis. Toute norme est une généralisation de la notion de 
valeur absolue d’un nombre et doit satisfaire aux relations 


fu+vl<lel+fvl,  eul=ftaetiul, 


où & est un nombre. 

Nous n'avons supposé nulle part que les fonctions dis- 
crètes u,/f figurant dans (94), (98), (100) sont des scalaires. 
Si le problèmo initial consiste dans l'intégration d'un 
système d'équations différentielles, dont la solution est un 
système de fonctions, c'est-à-dire un vecteur fonction, il est 
évident qu'alors u, f, lu sont des vecteurs fonctions dis- 
crètes. On doit en tenir compte lors de la définition d'une 
norme. 
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Ainsi pour un problème quelconque linéaire aux diffé- 
rences du type (94) la convergence 


Iu—U]—-0 pour +, k—+0 (101) 
découle de l'approtimation 
HU —f[-+0 pour +, h—0 (102) 


et do la stabilité 


LEA RAF AE 


Cette dernière relation signifie que, pour / quelconque, on 
a l'estimation suivante pour la solution correspondante « 


[ul <const ||f |}, (103) 


où const ne dépend pas de + et . 

Il est évident que si le problème est stable, (101) et (102) 
tendent vers zéro de la même manière, c'est-à-dire la pré- 
cision do la solution et l’approximation sont du même ordre 
de grandeur. 

La vérification de l'approximation, même dans le cas de 
problèmes compliqués, ne présente pas de difficultés. Vu la 
régularité de la solution exacte et compte tenu des équa- 
tions différentielles auxquelles elle satisfait on peut calculer 
la grandeur {U — f tout comme dans l'exemple simple 
envisagé ci-dessus. 

Les particularités du problème jouent un rôle important 
dans l'étude de la stabilité. Il est rare que l'on puisse, comme 
dans l’oxemple étudié, estimer directement «# en fonction du 
deuxième membre f. Pour un grand nombre de problèmes 
il vaut mieux vérifior la stabilité à l’aide de solutions parti- 
culières, la méthode généralisant ayant permis dans l’exem- 
ple cité de faire apparaître l'instabilité. Nous reviendrons 
à cette question ci-dessous. 

Passons aux problèmes non linéaires. Pour se faire une 
idée des difficultés apparaissant ici, nous allons commencer 
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par le problème suivant 
oU au \2 
+) = F (x, ” 
U(x, 0) =®D(x). 


Pour trouver la solution il est tout naturel d'utiliser le sché- 
ma aux différences 


(104) 


nt un PA be Led ne - 
T h Li (105) 
0 
Ur = Pr. 


Ceci se base évidemment sur le fait que le problème (104) 
peut être approximé par le problème (105), la solution U du 
premier satisfaisant presque le second. En effet U (x, ?) étant 
lisse, en utilisant l'équation (104) on obtient 

Der Un Un (Ur 9 

+ (—) = A+ 0 D, À 406) 


UF = qn. 


T 


Pour estimer l'écart de & à U on substitue dans (105) 
u = U + ôu et compte tenu de (106) on obtient 


Gun +1 Ga +9 Uri Ur | ôuR 41 — Our 
T h h 


Guy sy —Ôu 12 
+ (4) = O(t,h), 


+ 
(107) 


Oui — 0, J 


La présence d’un terme quadratique ne permet pas d'estimer 
d’une manière efficace le comportement de ôu. Supposons 
cependant qu'il y ait convergence et que ôu soit petit, par 
exemple ôu & U. Dans ce cas le comportement de ôu sera 
déterminé par les termes linéaires principaux de (107), 
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c'est-à-dire par l'équation 
n+i n n n n n 
ôu, — ôu, po Ti TU Guy, —ôuy 
T h k 


— O(x, h) (108) 


qui, à un coefficient près, coïncide avec l'équation envisagée 
antérieurement (82). Pour que cette coïncidence soit complè- 
te et que l’on puisse utiliser les résultats déjà obtenus, nous 
nous limileront à l'étude de (108) au voisinage du point 
considéré x,, {" où l'on à 


— 
L 


Usa Ur aU \n 
2H tn 2 (Te) 
ha x Jh 
et (108) peut être remplacé par (82), soit: 
Gun +1 Ôu}, ôup y — Ou} 
T h 
Comme nous l'avons déjà établi, la solution de cette équa- 
tion ne reste une grandeur petite qu'à condition que (84) 
soit vérifiée, celle-ci signifiant pour le moment que 


—=0 (rx, h). 


LL 
go ns Cp (109) 
Ainsi, pour que l'hypothèse de la convergence ne donne pas 
de contradiction il faut que la condition (109) soit vérifiée. 

On ne peut affirmer que la question de la convergence pour 
le problème (105) soit complètement résolue, néanmoins la 
condition (109) donne beaucoup. En particulier, si la solu- 
tion exacte a, dans un certain domaine, une dérivée posi- 


Q oU ? e e ? A 
tivo—>0, les équations aux différences ne peuvent ètre 


iiliséos: Pour vérifier la condition (109) au lieu des valeurs 
U, qui sont inconnues, on doit de toute évidence prendre 
u}, obtenues dans le calcul. 

Appliquons les raisonnements faits au cas général 
d'un problème non linéaire aux différences que l'on peut 
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écrire 
L (uw) = 0. (110) 
L'opérateur Z étant non linéaire, il n'y a pas lieu d'isoler 
le deuxième membre f. 
Avant tout on vérifie l'approximation, c'est-à-dire que 
la condition 


L(U)—0 pour t,k—0 (111) 


soit vérifiée. Si cette condition n'ost pas vérifiée rien ne 
donne d'espérer qu'il y aura convergence. 

On aurait pu généraliser la notion de stabilité aux 
problèmes non linéaires, en définissant celle-ci à partir de 
la relation 


u—U=£(u) — £(U). (112) 


Cette dernière no pouvant être vérifiée, ce n’est pas la peine, 
bien que pour la convergence c’est justement la relation (112) 
qui manque; en effet, vu l’approximation, la différence 
£ (u) — £ (U) cest petite. 

L'opérateur £ est non linéaire et ses propriétés lorsqu'il 
est appliqué à différentes fonctions peuvent varier. Natu- 
rellement, ce qui nous intéresse en premier lieu, ce sont 
des fonctions voisines de la solution exacte. Si, ici déjà, 
les propriétés de Z ne sont pas satisfaisantes, il n'y a pas 
lieu do s'attendre à la convergence. Mais pour u — U = ôu 
petits on peut obtenir l'information nécessaire sur Z (u) 
en considérant sa partie linéaire principale, c'est-à-dire 


L (U) + £'(U) ôu, (113) 


où £’ (U) est un opérateur linéaire (variante de £) appliqué 
à Ôu et dépendant de U comme d'un paramètre. Remar- 
quons que nous ulilisons ici le fait que Ôôu est petit, mais 
du point de vue formel ça n’a rien à voir avec la petitesse 
de T, À. 

Ainsi, en substituant dans (110) u — U + ôu, en rem- 
plaçant Z (U + Ou) par (113), compte tenu de (111), on 
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obtient 
£'(U)ôu—+0 pour t,h—+ 0. (114) 


Si ce problème provenant de l'opérateur linéaire Z’ (U) 
est stable, on peut s'attendre à ce qu’il y a convergence. 
En cas d’'instabilité, la convergence n'a pas lieu. 

En considérant l'opérateur £” (U) localement, dans un 
pelit domaine du plan x, {, où U (x, t) varie peu, on peut 
encore simplifier le problème en le ramenant à l'étude de la 
stabilité d’un problème linéaire à coefficients constants. 
En simplifiant, en simulant le problème il est important 
de ne pas omettre aucun de ses traits essentiels. T1 faut 
uno certaine habileté pour tenir compte au maximum du 
caractère spécifique du problème étudié. 

Résumons donc. Lors de l'étude de la convergence d'un 
problème aux différences il y a lieu, avant tout, de vérifier 
les conditions d’approximation. En cas d'un résultat posi- 
tif on procède à la linéarisation du problème, puis l'on 
maintient les coefficients constants, ce qui nous ramène à 
l'analyse do la stabilité d’un schéma linéaire aux diffé- 
rences à coefficients constants. Brièvement cette simplifica- 
tion, cette simulation du problème peut être représentée 
par la formule suivante 


L (u) — L'(U) ôu —+ lôu. (115) 


Problèmes 


1. Démontrer, par estimation directe, que le probli- 
me (108) est stable lorsque la condition (109) est vérifiée. 
2. Trouver les conditions de stabilité (de convergence) 
pour les schémas aux différences approximant les problèmes 
oU aU - 
ar Us: U (0, x) =U(x) 
et 
IU 
=F(U),  U(0)= Un. 
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3. Elaborer les schémas aux différences approximant les 
problèmes 


oU Ge) oU 
aan), U(, += U(r) 
et 
QU Op (V 
FT ET =0, U(0, x)-=Uo(x), 


oV oU 
ta =0, U (0, x) = Vo(x). 


Effectuer la linéarisation des équations aux différences 
obtenues. 


$ 6. CRITÈRE SPECTRAL DE STABILITÉ 


Pour les problèmes linéaires aux différences du type 
lu — f la stabilité signifie que # = f c'est-à-dire quo la 
solution ot le second membre sont pour +, À —> O0 du mêmo 
ordre de grandeur. Nous allons nous limiter à l'étude des 
opérateurs ? à structure en couches du type 

— n=0, 1,..., 
Lu — 


T 


(116) 


uw, 


Ici #” désigno la fonction discrète sur n-ième couche, c'est-à- 
dire l’ensemble u% pour # donné, À étant un certain opéra- 
teur linéaire transposant une fonction sur une couche en 
une fonction sur une couche et dépendant des paramètres 
T, k. Pour les opérateurs ? du type (116) le problème !u = f 
pout s’écrire sous la forme 


uni Nu +rf", n—=0, 1,..., 
u = v, | 


(117) 
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où /", v sont des fonctions discrètes données sur les couches. 
C'est là la stratification de l'opérateur Là savoir les fonc- 
tions discrètes ut, u%, ... peuvent être obtenues successive- 
ment à l'aide du même opérateur R, c'est-à-dire à l'aide 
de l'opérateur de transition d'une couche à l'autre. 

J1 y à lieu de trouver la relation existant entre les gran- 
deurs uw et f, v. Pour quelconque les formules (117) per- 
mettent d'exprimer «" en fonction de f" (n = N — 1, 
N —2,...) et v. En effet, 


uù = tft RuN-! = 
ne +R + __. = 


_;f NL pi - L Ref" Rs SON Ne — 
N-1! 

=T ÿ Vi D EUNES 00 (118) 
m=0 


où À" désigne la m-ième puissance de l'opérateur R et est 
lui aussi un opérateur linéaire. 

Pour estimer les grandeurs des fonctions discrètes sur la 
couche v, uw", f", Rv, ..., introduisons une norme quel- 
conque (par exemple {|| &" || — max |u}}). Comme uw" 

(1 


s'exprime en fonction de /", et v au moyen de 7", il est 
évident que la relation existant entre les grandeurs uw 
el f", v dépendra des propriétés métriques des opérateurs 
R", autrement dit de la variation de la norme des fonctions 
discrètes lorsque l’on applique AR”. Supposons que les 
opérateurs À" soient tels que, pour une fonction discrète 
quelconque sur la couche v, on ait l'inégalité 


RL I pm Ho | (119) 


où ph sont des nombres qui ne peuvent être diminués sans 
infraction à (119) ne serait-ce que pour une fonction v (p,, est 
appelé norme de l'opérateur R°"). 
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stimons la grandour #N en utilisant (118), (119) 


N-1 
pr D RE "14 Ro 


N-! 
LT 2 Pi + px Île 
ni 
<iN max Dm max 14 ex tel (120) 


Comme seules nous intéressent les valeurs finies de f{, on a 
0OSmr<N<t & OLm<N LIT. 


Il y aura stabilité si les coefficients auprès de max || f” || 
et ||[v || dans (120) resterons bornés pour t, À —+ 0. Comme 
TN <t, la condition de stabilité s'écrira 


Pm const pour T, À +0, mt < 1. (121) 


Il est évident que const ne doit pas dépendre de +, À bien 
que ph tout comme 2? et Z?" dépendent de ces paramètres. 
C'est là l'essence même du problème car pour tous t, À, m 
donnés la grandeur p,, est naturellement finie. 

Ainsi, pour l'opérateur ? de la structure stratifiée (116) 
la question de la s{abilité du problème aux différences revient 
à l'estimation des normes des puissances de l'opérateur R, 
‘c'est-à-dire à la vérification de la condition (121). 

Le problème envisagé au $ 4 donne, de toute évidence, un 
exemple de l'opérateur Z de la structure stratifiée. Dans ce 
problème nous avons d’un côté établi par estimation directe 
que (121) est vérifié sous la condition —1 << at/h < 0 et 
d'un autre côté, nous avons fail apparaître l'instabilité 
à l’aide de la solution particulière u, — (—1)*. Dans les 
cas plus compliqués il est rare que l’on arrive à vérifier (121) 
alors que l'on peut faire apparaître l'instabilité par générali- 
sation à un grand nombre de problèmes. 
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Pour étudier le comportement do R”"ven fonction de m”, 
le plus simple est d'étudier les fonctions propres de l'opéra- 
teur Z?, c'est-à-dire des fonctions discrètes sur la couche v 
qui, lorsqu'on leur applique Z?, se trouvent multipliées par 
des nombres À (appelés valeurs propres) 


Ro = Av. (122) 


Pour les fonctions propres 2”v = À"v, lorsque leur nombre 
est suffisamment grand, la grandeur À” permet de caracté- 
riser la norme p,. 

Considérons les fonctions discrètes sur la couche v — 
— {v,} déterminées et bornées pour toutes les valeurs de 
l'argument discret À, —o0o << k << oo. Supposons que l'opé- 
rateur linéaire Z? soit sur ces fonctions donné par la formule 


(Ro)u = D Gp Vntp k=0, +1, +2,..., (123) 
P 


où æ&h sont des coefficients donnés dépendant des paramè- 
tres t, À; p prend un certain ensemble de valeurs. En parti- 
culier, l'exemple envisagé au $ 4 est obtenu pour &, — 1 + 


at — AT 
+ +, & = —7—, pour les autres p on à a, = 0. 


Les fonctions propres de l'opérateur À du type (123) 
doivent satisfaire à la condition (122), c'est-à-dire 


D ŒpUh+p — AU k=0, Æ 1, = 2, …. (124) 
p 
Nous allons chercher la solution de cette équation linéaire 
aux différences sous la forme suivante 
Un = Vols (125) 


où g est un certain nombre et v, un facteur de normalisation. 
En substituant (125) dans (124) on obtient après simplifica- 
tion par og 


2 leg? = À, (426) 
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c'est-à-dire que pour g quelconque la fonction v (125) satis- 
fait à (122), avec À = À (q) (126). De ce grand nombre de fonc- 
tions nous choisirons seulement les fonctions discrètes (125) 
bornées (par rapport à X).Si]q| > 1, on a | v, | — oc pour 
k —+ oo ou pour À —> —oco. l’ar conséquent |gq | = 1. 
Quant à l'équation aux différences proprement dile nous 
l'envisagerons seulement dans le domaine réel. Cependant 
toute fonction réelle peut s'écrire comme une combinaison 
de fonctions complexes. C'est pourquoi l'utilisation de ces 
dernières peut donner une information utile sur les proprié- 
tés métriques de l'opérateur Z? dans l’espace réel, en parti- 
culier, sur sa norme. Dans le cas présent nous avons en 
tout deux fonctions propres, pour qg = 1 et q — —1. Le nom- 
bre de fonctions propres complexes est bien plus important 
el les propriétés de l'opérateur apparaissent sur ces fonc- 
tions. Posant qg = e' on peut écrire (125), (126) comme suit 


Un = Voeix?, (127) 
À = 3; apeirr, (128) 
p 


Ainsi, à chaque @ sur l'intervalle (0, 2n) correspond 
une fonction propre v4 (127) de valeur propre À (128). Il est 
évident que la grandeur v, ne joue aucun rôle, on peut donc 
poser w, = 1. 

Comme pour les fonctions propres R"v — }”v, on a 


Ro | = || "0 I<max [A |" lv, (129) 


la grandeur max | À |" est pour ainsi dire la norme de 
l'opérateur À” sur lé système de fonctions propres. Ce systè- 
me est uno partie de tout l'ensemble de fonctions discrètes, 
c'est pourquoi la norme de l'opérateur 2?” ne peut être que 
supérieure à max |A |”, 


max |A |" <pe (130) 
ç 
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En comparant cette inégalité avec la condition de sta- 
bilité (121), nous arrivons à la conclusion que pour vérifier 
cette dernière il faut que 


max {A/[”"<const pour Tt, h—0, mt<&t. (131) 
w 


Les valeurs propres À, tout comme @&, en fonction desquel- 
les elles s'expriment par (128), dépendent des paramètres 
T, h. Comme m = 1/t — oo, la condition (131) équivaut à 


max|Al&1+0O(t) pour +, h +0. (132) 
1 


Dans le cas contraire |A |" — oo. Par contre, si | À | — 
— { Let, on a 


[A |" (1 Let)" _ et, 


Cette dernière grandeur, bien qu'elle soit finie, peut être 
assez grande. C'est pourquoi au lieu de (132) on envisage 
parfois une condition plus forte 


max |A|<1 pour +t, À—+0,. (133) 
) 


Nous avons ainsi obtenu la condition nécessaire de stabi- 
lité du problème (117) avec l'opérateur Æ? du type (123). 
L'ensemble des valeurs propres À est appelé spectre de l'opé- 
rateur et la grandeur max | À |, rayon spectral. Aïnsi la 
condition (132) ou (133) est appelée critère spectral de sta- 
bilité. 

Pour l'exemple du $ 4 l'égalité (128) donne 

T T_: 
A=Î+a-—a ei, 
c'est-à-dire le spectre est un cercle (dans le plan complexe À) 
de rayon | at/h | et de centre au point 1 + at/. Il est facile 
de se convaincre que | À | < Î{ seulement dans le cas où les 
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mêmes inégalités 
—1<a-<0 


sont vérifiées. 

L'avantage du critère spectral do stabilité est qu'il peut 
faciloment être généralisé à un grand nombre de problèmes 
plus compliqués, en particulier aux systèmes d'équations. 
Revenons au problème (116), (117) et supposons quo les 
fonctions discrètes w”, f", v soient des fonctions vectorielles, 
c'est-à-dire qu'elles soient déterminées en chaque point 
de calcul par plusicours grandeurs. Tous les raisonnements 
ci-dessus restent vrais et seule la conclusion doit être changée. 
A savoir, (123) étant maintenant une égalité vectorielle, 
&y, Sont des matrices carrées du même ordre que les vecteurs 
v,. En substituant 


Un = Vo} — voeih, (134) 
où v, est un vecteur dans (124), on obtient 


A 
» ae PUy == À Vo 
P 


qui est un systèmo d'équations linéaires homogènes par 
rapport aux composantes du vecteur v,. Pour que la solu- 
rion existe, il faut que le déterminant du système soit nul: 


[D ape —A1| 0, (435) 
P 


où Z est une matrice unité. La différence d'avec le cas scalaire 
est que maintenant le spectre se compose de plusieurs bran- 
ches À,, À, . . ., définies comme les racines de l'équation 
(135) au lieu de (128). 

Ainsi le critère spectral permet de réduire l'étude de la 
stabilité au problème de l'estimation de la grandeur des 
racines d'une équation algébrique. 

Cette équation donne dans le cas général seulement la 
condition nécessaire de stabilité. Cependant si pour une 
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certaine classe de fonctions, pour une certaine norme, le 
système de fonctions propres est complet, c'est-à-dire que 
toute fonction de cette classe peut être approximée par une 
combinaison de fonctions propres, le critère spectral est 
alors la condition suffisante de stabilité. Considérons par 
exemple le cas scalaire où l'opérateur Z? est du type (123), 
de plus nous nous bornerons aux fonctions discrètes w — 
— {u,} pour lesquelles la série 


Di une" = 1 (4) [(136) 


est convergente. Chaque fonction discrète uw, donne ainsi 
naissance à une fonction périodique w (œ) pour laquelle la 
série (136) est une série de Fourier, u, étant les coefficients 
de cette série. Ces derniers, comme on sait, sont donnés par 
les relations 


2x 
Up = _ w (@) efk® dp (137) 
et 
1 T 
Diuè=s | lw(p)lFdp. (138) 
k 0 


On peut directement se rendre compte de la validité des rela- 
tions (137), (138). A cet effet il y a lieu de multiplier (136) 


par ei"? dans le premier cas ot par w(p) — > u,,eim? 
m 


dans le second, et d'effectuer ensuite l'intégration. 
Appliquons à (137) l'opérateur Z? de (123), on obtient 


27 


(Aux —= >, Gp | WW (p) eik+P) € dp— 
p 0 


27 27 
— _ | iw (wp) e4? » GneiP® dp = _ | w (p) À (p) ef*? do, 
U p 0 
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ceci en vertu de (128). En comparant cette égalité avec (137), 
on voit que (/?u), sont les coefficients du développement en 
série de Fourier de la fonction &w (p) À (œ). Par conséquent 
pour ces cocfficients on a des relations du type (138), 
c'est-à-dire 


Lu (p) À (p) F d. 


n 


>, (Ru)i = _. 
k 


En mettant max | À (ç) |? en facteur devant le signe de 
p 
l'intégrale et en remplaçant l'intégrale restante par D «À 
k 


on obtient 
D (Ru) max |A (p) P 2 uk. (139) 
k y k 

Définissons la norme de la fonction discrète par l'égalité 


lul=Qun), (140) 


où le facteur À a été introduit pour que, à la limite, pour 
hk —+ 0, la norme conserve son sens. En utilisant (140), on 
peut écrire (139) sous la forme 


I Ru] max |A} [ul (141) 
P 


Cette dernière inégalité signifie que pour la classe men- 
tionnée de fonctions avec la norme (140), la grandeur max || 
n'est pas inférieure à la norme de l’opérateur À? ot par consé- 
quent le critère spectral est un critère suffisant de stabilité. 

Arrêtons-nous maintenant sur une question importante 
pour les applications pratiques du critère spectral de stabili- 
té. Dans tout problème réel de calcul le nombre de points 
de calcul est fini et l'indice Æ prend un ensemble fini de 
valeurs. Dans les points limites l'opérateur Z? ne peut 
conserver sa forme standard (123), ne serait-ce que vu l'ab- 
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sence d’un ensemble complet de grandeurs u,,,. En ces 
points l'opérateur 2? s'exprime à l’aide de formules spéciales, 
réalisant telle ou telle condition aux limites. Pour pouvoir 
se faire une idée à quel point cette déformation peut changer 
les propriétés de l'opérateur À (123) standard, « sans limi- 
tes», il y a lieu de revenir à l'exemple du $ 4. 

Supposons que Z? soit donné sur des fonctions discrètes 
un (x = 0, 1, 2,..., K) par les formules 


(Ru), = (14 a +) TT TRE TS k—0O, 1,..., K—1, 


h 
(Ru)x =(. 
(149) 


Il est évident que ceci correspond à un problème différentiel 
sur un intervalle fini 0 < x < X — Kh avec une condition 
aux limites nulle sur l'extrémité droite, dont la solution 
pour a << 0 existe et est unique. 

Nous allons trouver les fonctions propres de l'opérateur À 
(142). Posant Rv = Àv on arrive à un système d'équations 
linéaires homogènes par rapport à Up, Uy, + +. Un! 


(a—1—a +) oi +a vu = 0, k—0,1,...,K—1, 
Avr = 0, 


Il est facile de voir que ce système a deux solutions non tri- 


viales 
T h 
{ La — 
Fe 
T 
O — 


he 


Un = pour À —0 


et 


Vo=Â, U=vo—...=Vvx=0 pour A=t+a—. 


Bien que la seconde solution entraîne la nn 1 + a+ < 


< 1, il est clair qu'un nombre aussi restreint de fonctions 


UB 


propres ne permet pas de mettre en évidence les propriétés 
de l'opérateur. 

En même temps il est peu probable que la déformation 
de l'opérateur en un point limite puisse notablement changer 
ses propriétés, apparaissant d’une manière si évidente sur 
les fonctions e“®, Ces dernières ne satisfont pas à la condi- 
tion ux — 0, c'est pourquoi nous allons les corriger, c'est-à- 
dire poser 


Un = ei, k=—0, 1, .., K—1, 
) cu 


Ur = 0 


et leur appliquer l'opérateur À (142). On voit aisément que 
pour tous Îles Æ << À — 1 la fonction (143) se transformera 
en elle-même, se trouvant multipliée par le facteur 


A=ita—a—eir, (144) 


et ce n’est qu'au voisinage du point limite que cette loi 
cesse d’être vraie. Appliquons de nouveau l'opérateur À (142) 
à la fonction obtenue, c'est-à-dire étendons la correction 
également au point À — 2. En répétant ce processus on 
obtient pour la fonction v, (143) 


(R'v), = Av, pour k=0, 1,..., K—m. (145) 


Ce qui nous intéresse c'est R" pour Tt, À —+ 0, c’est-à-dire 
pour m — Â/t—+ 00, K =1/h—> 00. Si alors h/t resto 
borné, en d’autres termes si »7 tend vers l'infini pas plus 
rapidement que X, il y aura toujours un intervalle de va- 
leurs k, où (145) est vérifié. Par conséquent l'estimation 
inférieure de la norme || 2” || en fonction de ||”, où À est 
donné par (144) c'est-à-dire est une valeur propre de l'opé- 
ratour non pertubé, resto vraie. 

On peut raisonner autrement, en considérant au lieu de 


(143) la fonction 
= (1 — —) env, (146) 
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En utilisant (142), (144) on trouve que dans ce cas 
(Ro — Av), = eith+ 0, k—0, 1,..., K—1, 


(Ru — DPF = 0, 


Comme XX = X — const, les dernières égalités signifient 
que 
iv — Àv = 0 (x) (147) 


bien que la fonction v (146) soit continue (v, — 1). On peut 
dire que v (146) est une fonction « presque propre » et À (144). 
une valeur « presque propre » de l'opérateur R (142). Nous 
arrivons de nouveau à la conclusion que les conditions aux 
limites doivent être considérées comme une perturbation de 
l'opérateur, pour laquelle la plupart de ses propriétés se 
conservent. Il est évident qu'en se donnant des conditions 
aux limites absurdes on peut abîmer l'opérateur standard, 
mais on ne peut pas corriger ses défauts par des conditions 
aux limites appropriées. 

Dans le cas général l'étude de ces problèmes n'est pas 
simple. Cependant les raisonnements ci-dessus nous mon- 
trent que pour les problèmes sur un intervalle limité la 
condition nécessaire de stabilité reste l'application du 
critère spectral à un opérateur standard limité. 


Problèmes 


1. À partir de la définition de la norme de l'opérateur 
(119) on a 


Ru | = RAR Ke RTE. pre], 


c'est-à-dire p, < p". Par conséquent pour que (121) soit 
vérifié il suffit que p® < const, ce qui équivaut à la condi- 
tion suivante pour la norme de l'opérateur /? 


Pi <1+0 (1), (148) 
qui est par là même la condition suffisante de stabilité. 
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Montrer que pour l'opérateur 7? de la forme (123) la con- 


dition (148) pour la norme || u || = max |, | équivaut à la 
k 
condition 
D lapl<1+0(). (149) 
p 


Pour l'opérateur 
at 
(Run = Un — oi (uns —Un-1) 


étudier la stabilité à l’aide de (149) et du critère spec- 
tral. Comparer les résultats obtenus. Définir le problème 
différentiel correspondant à cet opérateur. 

2. Pour le problème 


U 0V 

+ = 0, U (0, x) =Uo(r), 
4 OU 

+0 V (0, x) = Vo(x) 


trouver les différents schémas aux différences ot étudier leur 
stabilité. 
3. Elaborer le schéma aux différences pour la solution 

du problème 

5184 ôf (U ô oU 

ET CCR (1) EE U (0, x) = Uo(x), 
où f(U), u(U) sont des fonctions données. Etudier 
l'approximation et la stabilité (sur le modèle linéaire). 


Ê TABLISSEMENT DES FORMULES 
$ 7. De cALcuL 


Maintenant que nous connaissons les conditions aux- 
quelles doit satisfaire le problème aux différences nous allons 
passer à son établissement. Dans chaque exemple concret 
nous avons simplement remplacé chaque dérivée entrant 
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dans une équation différentielle initiale par la relation aux 
différences correspondante. Mais ce n’est pas là la méthode 
la plus courte, ni la plus simple. 

La composition du problème aux différences commence 
par le choix du réseau de calcul, c'est-à-dire d'un ensemble 
discret de points remplaçant le domaine continu des varia- 
tions des variables indépendantes. 

En principe cet ensemble est arbitraire. On peut augmen- 
ter la concentration des points des parties les plus importan- 
tes afin d'y obtenir une précision plus grande. On peut se 
donner la loi de formation du réseau en fonction de la solu- 
tion obtenue durant le calcul. Mais sauf prescription spéciale, 
il est préférable de prendre le réseau régulier, déterminé 
par un nombre minimal de paramètres. Ceci facilite notable- 
mont l'étude du problème aux différences. 

Supposons que le réseau de calcul ou la loi de sa forma- 
tion soient donnés. Si l'irrégularité de réseau est peu accu- 
sée, c'est-à-dire ses paramètres varient peu d'un point à 
l'autre, sur de petits domaines il peut être simulé d'une 
manière uniforme. Ultérieurement nous considérerons un 
réseau régulier (pas obligatoirement rectangulaire) déterminé 
pour les problèmes à deux variables indépendantes x, t{ 
par deux paramètres seulement, à savoir les pas X,T du 
réseau. Sur ce réseau on détermine les fonctions (ou les 
vecteurs fonctions) u%, f#, @n, . .. qui sont les fonctions 
discrètes des arguments discrets kÆ, n, (numéros des points). 

L'étape suivante est l'établissement des équations aux 
différences, c'est-à-dire des relations arithmétiques entre les 
grandeurs +, k, ul, ff, Œn, . . . Comme nous nous basons 
sur le principe de la convergence de la solution du problème 
discret vers la solution de l'équation différentielle, les 
équations aux différences doivent satisfaire à certaines 
conditions. Ci-dessus nous les avons formulées comme des 
conditions d'approximation et de stabilité. 

En utilisant les formules de la dérivation numérique en 
vue de remplacer les dérivées par des différences finies, il 
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est facile d'écrire telles ou telles relations qui à la limite, 
pour une fonction lisse quelconque, deviendront les équa- 
tions différentielles initiales. Mais ici le succès n'est pas 
immédiat car la majorité des schémas logiquement possibles 
sont instables. 

La recherche d'un schéma aux différences convenable 
peut être plus cfficace si l'on utilise un artifice que nous 
ex poserons tout d'abord sur le même excmple simple 


au ou 
a Te ge 0 | 450) 


U (x, 0) =U, (x). 


Nous allons fixer la forme de la naille de calcul, c'est-à- 
dire indiquer les points où les valeurs de la fonction discrète 
sont celles que nous voulons lier aux relations aux diffé- 
rences. Pour le problème étudié prenons sur le réseau x; — 
= kh, "= n7t quatre points dont les numéros sont (k, x + 1), 
(4 — 1, n), (k,n), (k + 1,n) formant une telle maille 
(fig. 7). 

Lo problème initial (150) est linéaire. Il est donc tout 
naturel de chercher les équations aux différences sous la for- 
me de relations linéaires. La forme générale d’une telle 
relation entre les valeurs de la fonction discrète aux points 
mentionnés est 


{ 
unti= our, + our + œur,, — À, au, (151) 


Posant uf = U, (x,) on obtient pour tout ensemble @&; un 
certain problème aux différences. Nous allons choisir les 
coefficients &, tels que le problème (151) donne une approxi- 
mation du problème initial et soit stable. 

Si pour vérifier la stabilité on a l'intention d'utiliser 
le critère spectral, il y a lieu d'écrire le problème (151) 
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nel 


Uk 
T 
h h 
HE u£ Up 
Fig. 7 
sous la forme usuelle !4 = f en posant 
un +1 —ÿ, Apple 4 p 0 
p Ù 
D ———_— = 15 
de . Ph Us (tn). GE. 


uÿ, 


On sait que la condition nécessaire de stabilité s'écrit dans 
le cas présent sous la forme de l'inégalité suivante 


[5 œpeirr|<1. (153) 
P 


Nous avons ainsi une condition imposée aux coefficients &;. 

Passons à l'approximation signifiant que !U — f +0 
pour Tt, À — 0. Pour trouver l'ordre de grandeur de IU — }, 
utilisons comme toujours le développement de. la solution 
de (150), à savoir de la fonction U (x, ft), en série suivant 
les puissances de tv, À, au voisinage du point central de la 
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cellule de calcul x,, {”. On a 


sé 2h2 33 
Urip = U + phU, + je Uxx + U rex + 


n 2 3 
Ur =U ++ Ur + Ur + .., 
(154) 


où U, Uy, ... sont pris au point central. Comme U (x, t) 
est la solution de l'équation (150) pour cette fonction on a les 
égalités 
U, = — aU,, 
Ur = U x, 
Ur = — QU xxx, 


(155) 


obtenues par dérivation de (150). 
En utilisant (154), (155) on peut trouver la combinaison 
qui nous intéresse des valeurs de U: 


UE — S apUr4p— (1 + Gp) U — (ra+h > Pan) Ux+ 
P 


++ (rate 5 Pa) Urx— 


+ (ra + hé » pop) ÜUxxx + (156) 


ne différant de LU — f que par le facteur 1/+. L'ordre de l’ap- 
proximation dépendra des premiers termes non nuls de ce 
développement. Les grandeurs U, U,, U,,, ... doivent 
être supposées quelconques et indépendantes. En annulant 
les coefficients de ces grandeurs, on obtient une chaîne 
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d'équations 


ta+h D pas =0, 
ra RS pa, =0, (157) 


Ta +R D pan = 0, 


qui sont les équations permettant de déterminer æ&, dont 
l'ordre, par rapport à t,:, augmente. Plus grand est lo 
nombre d'équations auquel on peut satisfaire, plus l'ordre 
do l'approximatlion sera élevé. 

Nous n'avons à notre disposition que trois coefficients 
non déterminés &,. C'est pourquoi nous pouvons satisfaire 
à trois des équations (157) au plus. En désignant at/h par 
r, on peut écrire ces équations comme suit 


Li + a+ —=1, Luy—d—=Fr, Qu +a= rt, 


En les résolvant, on obtient 


L_y — Dr. Xo—=1—7r?, a = TE . (158) 
Pour ces valeurs de «, l'ordre de grandeur du premier terme 
différent de zéro dans le développement (156) est 1°, rhÿ. 
Par conséquent ZU — f = 0 (x?, h°?), c'est-à-dire l'approxi- 
mation est du second ordre par rapport à T, . 

Il ne reste qu'à satisfaire la condition de stabilité (153). 
En substituant (158) dans (153), on obtient : 


2 2 _— . 

D) a peire = — 2 EPA — TE OT pig — 
= 1—r2+r2cos p— ir sin o — 
CE, .. P _P 
= 1 — 2r2sin 2 i2rsin 7 COS + 
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Par conséquent 


| D, apeirg = (1 — 2r? sin? 2)" + (2r sin Ÿ cos &)" = 


— 14 — 4r2sin2 L + 47 siné L + 4r2sin2 -L cos P — 
2 2 2 2 
— 1 —4r2 sint À y + 4r* sin + == 


me À — 472 (1 — 72) sin“ +. 
Il est évident que pour vérifier l'inégalité (153) il faut que 
r L'1, c'est-à-dire 


ja <1. (159) 


C'est la condition de stabilité du problème aux différences 
(151), (158). 

Si pour trouver les coefficients &, on se borne aux deux 
premières équations (157), le schéma obtenu sera de toute 
évidence une approximation du premier ordre © (t, k). 
Comme ces deux équations contiennent trois inconnues &, 
il y a un grand nombro de schémas pouvant être utilisés. 
L'un de ces schémas a été envisagé au $ 4. Si l’on utilise 
seulement la première des équations (157), l'approximation 
ne se trouve pas assurée, car dans ce cas l'erreur sera finie. 

Il va de soi que la méthode décrite peut également être 
utilisée pour trouver les formules de calcul d’autres problè- 
mes. Avant tout, un raisonnement à priori nous permet de 
choisir la forme de la maille de calcul et celle des relations 
aux différences contenant dos coefficients indéterminés. En 
exigeant ensuite que soient remplies les conditions d'appro- 
ximation et de stabilité, le problème revient à la solution 
d'un système d'équations et d'inégalités algébriques. 

Lorsque l’on utilise la méthode des coefficients indéter- 
minés pour la solution des problèmes compliqués on est 
obligé de faire des calculs analytiques compliqués. 
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Pour simplifier les calculs il y a lieu de renoncer à la 
généralité et de se borner à la solution des formes simples 
d'équations aux différences. Par exemple dans le problème 
envisagé, on aurait pu au lieu de (151) chercher le schéma 
aux différences sous la forme 


un+1 = uy 


6, ee LC uy — _ , (160) 


T 


les coefficients c,, c, étant indéterminés. 

Il est évident que si les équations aux différences satis- 
faisant aux conditions imposées (forme de la maille de 
calcul, ordre de l’approximation, forme de la relation aux 
différences) existent, elles peuvent toutes être obtenues par 
la manière décrite. 

Nous allons aborder encore une question. Lors de l'étude 
du schéma aux différences (151) nous l'avons écrit sous la 
forme lu = f, de plus, nous avons déterminé ! par la for- 
mule (152) et pour une raison « mystérieuse » divisé (151) 
par +. Il semble plus naturel de définir ! par l'égalité suivante 


lu = ut! — ù AP pe (161) 


p=-i 


Remarquons que dans ce cas (pour les mêmes «,) on obtien- 
dra LU — f = 10 (x°, hk?) c'est-à-dire une approximation 
d'ordre plus élevé. Cependant les propriétés du problème 
ne dépendent pas, de toute évidence, des désignations, des 
notations et de la méthode d'étude. L'avantage d’une meil- 
leure approximation est illusoire et se perd vu la stabilité 
plus faible. En effet, comme il est facile de le voir, pour 
le problème !u = f où ! est donné par (161) on obtient 
l'estimation uw — f/x, ce qui, conformément à notre défini- 
tion, signifie qu'il y a instabilité. Néanmoins ceci n’a pas 
d'influence sur la convergence car on a un ordre de réserve 
dans l’approximation. 


78 


La remarque faite montre la non-unicité formelle appa- 
raissant lorsque l'on divise le problème de la convergence 
en approximation et stabilité. 

Examinons encore une méthode permettant de trouver 
les formules de calcul. Cette méthode a un domaine d'ap- 
plication plus restreint mais se trouve être souvent suffi- 
samment efficace. 

On peut dire que notre méthode d’approximation des 
problèmes différentiels est locale. Lors do l'établissement des 
équations aux différences toute notre attention se porte sur 
une maille de calcul, sur le domaine voisin du point de calcul 
dont les dimensions sont de l'ordre de +t, k. Or, de même que 
toute fonction lisse peut localement être supposée linéaire, 
tout problème envisagé dans un petit domaine où la solu- 
tion varie peu peut être approximé par un problème linéaire 
à cocfficients constants. 

Par exemple au lieu de l'équation 


ou | aU 
on 0) 0, Ge 
on peut au voisinage du point +, {" considérer l'équation 
oU oU 
a + 0 (163) 


où aë — a (xx, !", Uf). En effet, si l'on substitue la solu- 
tion (lisse) de la première équation dans la seconde, celle-ci 
est satisfaite à © (xt, X) près car 


(a(x,t, U(x, t))— a?) + = O (x, h). 


Ajoutons à ce qui vient d’être dit que l'étude de Ja sta- 
bilité des problèmes aux différences ($ 5) que nous avons 
faite ci-dessus est basée sur leurs modèles linéaires. Toute 
la théorie de la stabilité concerne en fait la stabilité des 
équations linéaires aux différences. 
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Les raisonnements ci-dessus montrent qu'il y a lieu 
d'onvisager les méthodes d'établissement des formules de 
calcul en nous limitant aux problèmes linéaires à coeffi- 
cients constants. 

Ci-dessus, lorsque nous avons exposé différentes proprié- 
tés sur l'exemple du problème (150), nous avons systémati- 
quement évité d'utiliser le fait que nous connaissons la 
solution exacte et que celle-ci est donnée par une formule 


explicite, à savoir 
U (zx,1) = U, (x — at). (164) 


En l'utilisant nous risquions d'obtenir des affirmations 
vraies seulement pour les problèmes dont on connaît la 
solution, notre étude n'ayant plus alors de sens. Nous avons 
fait appel au problème (150) qui représente une certaine 
classe de problèmes et nous avons utilisé celles de ses pro- 
priétés qui sont typiques pour cette classe. Par contre 
l'existence d'une formule donnant l'expression explicite de 
la solution caractérise les équations linéaires à coefficients 
constants. 

Comme à l'heure actuelle ce sont justement des problèmes 
de ce genre qui nous intéressent, on est en droit d'utiliser 
cette propriété. Voyons la manière de procéder d’abord sur 
l'exemple de l'équation (163). 

Nous avons besoin d'obtenir une relation exprimant 
ut! en fonction des grandeurs de la »-ième couche up, 
uñ+1, . . . La formule de la solution exacte (164) appliquée 
au cas envisagé donne 


UÙ (rn, ET) = U (ax — ant, 2). (165) 


Par conséquent, pour trouver wf*} il faut connaître la valeur 
de la solution au point x = x, — a%t de la n-ième couche 
(fig. 8). Mais sur cette couche nous avons un ensemble discret 
de valeurs u}, uhzs, . . . c’est-à-dire la fonction discrète. 
Pour calculer la valeur dont nous avons besoin, procédons 
à l'interpolation. 
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Fig. 8 


Supposons que 
The LTh — TL Th (166) 


L'interpolation linéaire donne alors 


T T 
U" (Tr — Qt) = a? unit ( 1 — a +) un, (167) 
el conformément à la formule (165) on obtient 


n+i nn __ n T n _mn 
unti= ur — a} (u un _4), 


h k 


qui est l'équation aux différences que nous connaissons déjà 
et qui est stable pour 


T 
0<ar -—<1. 
On voit ainsi que la condition de stabilité dans le cas 
présent coïncide avec la condition (166) assurant l'interpola- 


bilité de la formule (167). Ce n’est pas étonnant car pour assu- 
rer la stabilité il faut prendre en considération le domaine 
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de dépendance de la solution. C’est là justement ce qui est 
essentiel pour la méthode proposée. 

Maintenant que les traits essentiels de la méthode sont 
clairs, sans aspirer à la généralité, nous allons les exposer 
pour le cas où lo problème initial consiste à intégrer un sys- 
tème d'équations différentielles aux dérivées partielles du 


type 
T9 (0) (168) 


pour des conditions initiales données. U est ici le vecteur 
fonction de x, t cherché et % un opérateur différentiel (par 
rapport à x) pour lequel le problème posé est correct, c'est-à- 
dire sa solution existe, elle dépend d'une manière unique et 
continue des données initiales. Ces problèmes sont dits 
d'évolution car ils décrivent le développement dans le temps 
d'un certain état initial. 

La première étape consiste à créer un modèle linéaire 
du système donnant une approximation locale du systè- 
me (168). L'opérateur Z est une certaine fonction vecteur 
donnée des arguments U, U,, U;,x, ... 


D(U)=f{U, Ux, Urx, ...). 
Donc au voisinage des valeurs UF, (U,)h (Uxx)# , - .. 


on peut approximer l'opérateur Z par l'expression linéaire 
suivante 


D (0) — 5 + Pur (U — UF) + (fu) (U x — (Ux)7) +... (169) 


fu, fu, - . . désignent les matrices des dérivées des compo- 
santes du vecteur f par rapport aux composantes des vecteurs 
U, U,, ...et fñ, U?, (U.)h sont les valeurs des grandeurs 
correspondantes au point 2, éd”. 

Supposant que la solution U (x, t) soit lisse, c'est-à-dire 
que les différences U — U?, U, — (U,)*, . . . soient petites, 
au lieu de (168) considérons le système d'équations différen- 


82 


tielles linéaires à coefficients constants 
2 DUC. (170) 


où D est un opérateur différentiel linéaire, s'exprimant par 
la partie homogène (169), et C une constante tenant compte 
des termes restant de (169): 


DU = (fu)U + (fu) Ux+ ..., 
C=f}— (ur Ur —(u,)} (0)? + 


La théorie des équations différentielles aux dérivées par- 
tielles nous apprend que la solution d’un système d'équations 
linéaires à coefficients constants (170) peut s'exprimer par la 
formule suivante 


U (x, t) = QU (x, 0) + Ct, (171) 
où Q est un opérateur intégral linéaire, à savoir 
QU'(x, = | QE DU(R+HE, O)dE, (172) 


et q la matrice des fonctions correspondant au système (170). 
Nous n'allons pas nous arrêter ici sur les méthodes permet- 
tant de l'obtenir. Mentionnons par exemple que pour l'équa- 
tion de la conductibilité thermique 


ou dU 


Wa,  a>0, (473) 
la fonction g est de la forme 
+2 
q(x, t)= a (174) 


2V nat 


Pour résoudre le problème (150) on peut également écrire la 
formule (164) sous la forme (171), (172), ceci à condition de 
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poser 
q(x,t) = Ô (x + at) 


où Ô (x) est la fonction delta. 

La formule (171) permet d'exprimer la solution à l'instant 
ati — fx à l'aide de la fonction U (x, {"). Dans notre 
cas cette dernière est représentée par la fonction discrète ux. 
Donc, pour utiliser la formule (171) trouvons tout d'abord la 
fonction d'interpolation Pu" (x). Comme toute notre étude 
porte sur le voisinage du point x,, {”*, nous écrirons cette 
fonction comme suit 

Pu" (x) = D ur, nPm (t), (175) 
où P,, (x) sont les polynômes correspondants. 

La fonction d'interpolation (175) peut également être uti- 
lisée pour la détermination des valeurs (U,)%, . . ., entrant 
dans l'expression (169), et par conséquent, dans les cocffi- 
cients D et le deuxième terme .du deuxième membre C 
de (170). 

En substituant Pu" (x) dans (172) on obtient 


QPu" (x) = 2 CAT ENS 


O0 


am= | QE, ©) Pm (tn +8) db, 


et en vertu de (171) on peut écrire 
url e D anthim + CT, 
m 
c'est-à-dire la formule de calcul dont on a besoin. Il est évi- 


dent que cette dernière est simplement une formule de 
quadrature pour (171), (172). 
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Nous allons faire plusieurs remarques. Le système (170) 
donne une approximation du système initial (168) à des 
termes du second ordre par rapport à Tt, À près, car l'erreur 
de l'approximation linéaire de (169) est du second ordre par 
rapport à U — Uf, U; — (U.,}$, ... Afin de simplifier 
l'opérateur linéaire D, on peut en écrivant  (U) sous la 
forme (169), ne tenir compte que des dérivées d'ordre élevé 
(comme dans (162), (163)) car la stabilité dépend essentielle- 
mont de la formo des lermes des équations aux différences 
correspondant aux dérivées d'ordre supérieur. Il est vrai 
que ceci a pour effet d'abaisser l'ordre de l'approximation 
jusqu'à © (t, h). 

Pour les mêmes raisons, lors du calcul des grandeurs U}, 
(U,)#, ..., entrant dans les coefficients et dans le second 
membre de (170), il n’est pas obligatoire d'utiliser la fonc- 
tion d'interpolation (175). A cet effet on peut utiliser des 
expressions quelconques aux différences donnant évidemment 
une approximation de ces grandeurs. 

Tous les schémas aux différences obtenus par la méthode 
exposée ne diffèrent les uns des autres que par la méthode 
d'approximation locale (169) et par la forme de l'interpola- 
tion (175). Si l'interpolation utilisée est suffisamment précise 
et permet de tenir compte du domaine de dépendance de la 
solution d'une manière efficace, le schéma aux différences 
est satisfaisant. Son étude s'effectue alors par les méthodes 
habituelles. 


Problèmes 


4. Quelle est l’imprécision, en tant qu’ordre par rapport 
à t, k, que l'on peut admettre lors de la résolution du systè- 
me (157)? . 

2. Trouver tous les schémas stables du type (151) et (160) 
pour le problème (150) du premier ordre par rapport à +, À. 

3. En utilisant la méthode des coefficients indéterminés 
trouver le schéma aux différences du type (151) pour la 
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solution de l'équation do la conductibilité thermique (173). 
Pour quelle relation entre + et À ce schéma donnera-t-il 
une précision maximale? 

4. En utilisant la formule de la solution exacte (171), 
(172),(174) de l'équation de laconductibilité thermique (173), 
construire le schéma aux différences, en utilisant pour 

= {": 

a) l'interpolation linéaire sur doux points 2,1, ÆZr+1: 

b) l'interpolation linéaire par morceaux sur %y-1, Zk 
pour æ <Z Zn Ot SUT ZT, Tr+1 POUT T> Ty; 

c) l’interpolation quadratique sur trois points zx-1 
Th: Th+1: 

Etudier l'approximation et la stabilité de chacun des 
schémas obtenus. 

5. La solution du système d'équations 


oU OV 
ôt EE Ôx = 0, 
9V oU 
HE À ôx 0 


est donnée par les formules 

U(x, t)=(U(x+t, 0)+U(x—t, 0)—V(x+4, 0)+ 
+V(x—t, 0)), 

Vx, D=+(V (a+, 0)+V (x—t, 0) —U (xt, 0)+ 
+U(z—t, 0)). 


Les utiliser pour construire le schéma aux différences cor- 
respondant au système 


20 0p09 56 


ôt ôzx ; 
1] 4 oU 
8 = 8 0 


où p (V) est une fonction donnée. 
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$ G. SCHÉMAS AUX DIFFÉRENCES IMPLICITES 


Lors de la résolution d’un problème concret par une 
méthode numérique il y a lieu de choisir les valeurs des 
paramètres de la méthode, c'est-à-dire les pas du réseau +, k, 
soit la quantité et la disposition des points de calcul. La 
solution de cette question dépend d'un grand nombre do 
facteurs différents : des propriétés de la solution, desexigences 
à la précision des calculs, des moyens disponibles de réalisa- 
tion de la méthode, c'est-à-dire de la puissance des ordina- 
teurs utilisés, etc. C'est pourquoi on ne peut donner une 
méthode pouvant être utilisée pour un grand nombre de 
problèmes. 

L'étude théorique de la convergence de la méthode n'est 
effectuée que pour résoudre des questions de principe. Pour 
les problèmes réels + et À ne peuvent tendre vers zéro. C'est 
pourquoi en choisissant le réseau de calcul de telle façon que 
la précision soit suffisante, on se base, en général, sur une 
certaine information sur les propriétés caractéristiques de la 
solution et l’on exige que le réseau soit suffisamment dense 
pour que ces propriétés puissent apparaître ainsi que les lois 
qui nous intéressent du comportement de la solution. A cet 
effet, en règle générale, on n'a pas besoin d’avoir un grand 
nombre de points. 

Cependant ce raisonnement n'est pas vrai pour toutes les 
méthodes numériques car il ne tient pas compte des proprié- 
tés individuelles, internes de la méthode. Sans toucher à tous 
les aspects de cette”question, nous nous arrêterons sur le 
plus simple. 

Toutes les équations aux différences concrètes que nous 
avons étudiées ci-dessus se sont trouvées être stables pour 
certaines conditions du type 


T T 
— <const, -7< const (176) 
imposées aux pas du réseau. 
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Jl est évident que lorsque l'on choisit les paramètres du 
réseau de calcul on doit en tenir compte. Souvent pour les 
valeurs de t, À donnant la précision requise les conditions 
(176) no sont pas vérifiéos. C’est pourquoi on est obligé de 
diminuer t (et ceci notablement), ce qui augmente consi- 
dérablement la quantité de calculs. 

Les conditions (176) étant exprimées à l'aide des para- 
mètres de la méthode numérique et non de ceux du problème 
initial ct comme elles caractérisent cette méthode, on peut 
essayer d'élaborer des méthodes imposant des conditions plus 
faibles aux pas du réseau ou mêmes sans ces conditions. 

Les schémas aux différences étudiés ci-dessus donnaient 
l'expression explicite en chaque point d’une couche tempo- 
relle donnée, u?*!, en fonction des valeurs de la solution en 
plusieurs points voisins de la couche précédente u%, uk41, 
c'est pourquoi ces schémas sont dits explicites. Pour ces 
schémas les conditions de stabilité du type (176) expriment 
qu'il est nécessaire de tenir correctement compte du domaine 
de dépendance de la solution. Cette relation n'est pas tou- 
jours aussi apparente que dans l'exemple du $ 4. 

Considérons le problème 


au o2U 
ts, U(, x) = Us (x), | 


— 00 Lr<oo, OLILT 


(177) 


et le schéma aux différences donnant son approximation 


un+i — u} up — Qu; +uÿ_, 
4 — = TR ’ u}, x U5 (ta), (178) 


k=0, +1, +2, ..., n=0,1,...,N. 


En utilisant le critère spectral, on peut facilement voir 
que la condition de stabilité du problème (178) s'écrit 


| 
T<T (479) 


Comparons les domaines de dépendance de la solution des 
problèmes (177) et (178). On sait que (voir (172) à (174)) 
pour l'équation de la conductibilité thermique (177) ce 
domaine est la droite —oo << x << oo toute entière. Mais 
dans le cas du problème aux différences (178) le domaine 
de dépendance pour les points de la N-ième couche seront 
les points de la couche zéro, remplissant l'intervalle de 
largeur finie 2WVk (fig. 9). Ainsi, bion que l’on ne tienne que 
partiellement compte du domaine de dépendance de la 
solution exacte, le problème (178) sous la condition (179) 
est stable. La contradiction se trouve évincée par le fait 
que pour T, À —+ 0, sous la condition (179), la largeur de 
l'intervalle mentionné croît indéfiniment. En effet, si 
t = T/N en vertu de (179) on a T/Nh? LL 1/2a et 


Nh> #7 —+ oo pour T, À —+ (0. 


Par conséquent à la limite on tient compte correctement du 
domaine de dépendance de la solution, ce qui de nouveau 
est lié aux limitations imposées par l'inégalité (179) aux 
pas du réseau. 

Ainsi, les schémas aux différences n’ayant pas de limita- 
tions importantes sur les pas du réseau, doivent être assez 
compliqués. Pour tenir compte d'une manière efficace du 
domaine de dépendance, il faut lors du calcul des grandeurs 
aux points d’une couche donnée utiliser un grand nombre 
do points de la couche précédente. Abordons la cause do 
l'instabilité et de l'apparition de la condition (179) d'un 
autre côté, quelque peu formel. 

En vérifiant la stabilité du schéma aux différences (178) 
à l’aide du critère spectral, nous étudions en fait le comporte- 
ment des solutions particulières du type 


un = Aeih®, (180) 
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Fig. 9 


c'est-à-dire utilisons la méthode de Fourier. En substituant 
(180) dans l'équation aux différences (178), on arrive à 


Anti àn ” _ MAS, (181) 


T 


où 
M =00 0%. (182) 


h? 


d'où l'on tire 
A A—iM)i = (A—TM)"i, (183) 
et comme seul le cas | À | < 1 nous intéresse, on obtient la 
condition | 1 — TM | < 1, c'est-à-dire (179). 
L’exposant nr dans (181) peut être interprété comme un 


indice (numéro du pas suivant t) et la relation (181) comme 
une équation aux différences correspondant à l'équation 
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différentielle 


= —MA,  A(0)=1, (184) 


où À est une constante positive, aussi grande que l’on veut, 
vu que À est petit. 
La solution exacto de l'équation (184) est 


A=e-Mt, (185) 


et la solution approchée donnée par (181) c'est-à-dire (183), 
convergera vers elle pour t —+ 0, car (1 — TtH4)t/ + e-Mi, 
Do ce point de vue la méthode (181) (c'est simplement la 
méthode d'Euler) d'intégration de l'équation (184) est 
acceptable. Mais cette méthode a un autre inconvénient. 
Alors que pour la solution exacte (185) l'estimation [A|<1 
est vraie pour M >> 0 quelconque, pour la solution approchée 
(183) elle ne l'est que pour des A7 potits satisfaisant à 
l'inégalité | 1 — TM |< 1. Mais c'est justement cette pro- 
priété de la solution qui nous intéresse. 

On peut facilement comprendre la cause de cet inconvé- 
nient : la solution exacte (185) est une fonction à décroissance 
rapide (pour M grand) et la méthode (181) utilise pour le 
calcul de la dérivée sa valeur disparaissant instantanément, 
non utilisable pour le pas suivant (—A/À" dans le second 
membre de (181)) (fig. 10). On peut rectifier la situation en 
prenant la valeur de la dérivée d'un pas à l'avance —M2"*1, 
c'est-à-dire 


Anti an 
T 


= — MM, (186) 
Dans ce cas au lieu de (183) on obtient 


An 1 
TR AA (87) 


Anti == 
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Ase M 


Fig. 10 


Cette solution, tout comme (183), converge de toute évidence 
vers la solution exacte pour T7 —>0, mais à la différence 
de (183) elle satisfait à la condition | À"*1 | < 1 pour tout 
M > 0 aussi grand que l’on veut, c'est-à-dire laisse bien 
voir la propriété essentielle de la solution exacte, à savoir 
sa décroissance rapide. 

Revenant au problème (177), on peut dire que la condi- 
tion de stabilité (179) est levée si on élabore un schéma aux 
différences correspondant à (186) et non à (181). En tout cas, 
la vérification de ce schéma sur des fonctions du type (180) 
ne fait pas ressortir d'instabilité. On a de toute évidence le 


92 


Fig. 11 


schéma suivant : 


TH —u} : nt} —Quftlpunti 
T du h2 } 
uk = Uo(xr), 559) 


k=0, +1, +2, ..., n=0,1,...,N, 


qui ne diffère de (178) qu'en ce que les grandeurs dans le 
second membre sont prises non pas sur la »#-ième, mais sur 
la (7 + 1)-ième couche; la maille de calcul correspondante 
est donnée fig. 11. I1 en découle que chacune des équations 
(188) reliant les valeurs u?*! en trois points ne donne pas une 
expression explicite pour w‘!, pour trouver ces dernières 
il faut résoudre un système d'un nombre peut-être pas 
infini, mais en tout cas grand, d'équations linéaires (188). 
Ces schémas sont dits implicites. Nous exposerons plus loin 
les méthodes de résolulion du système (188), maintenant 
nous allons continuer son étude. 

11 est évident que le problème (188), tout comme (178), 
donne une approximation du problème différentiel initial 
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(177). La vérification de la stabilité à l'aide de la fonc- 
tion (180) est basée sur le critère spectral. Nous l'avons 
formulé seulement pour des schémas explicites, son utilisa- 
tion pour des schémas implicites doit être argumentée. 
Cependant dans le cas du problème (188) on peut démontrer 
la stabilité directement. 

Modifions quelque peu le problème pour en faire le 
calcul moins ardu. Ainsi, nous allons considérer les équa- 
tions aux différences (188) seulement pour un ensemble fini 
(même très grand) de valeurs # et les compléter aux points 
extrômes par des conditions aux limites en donnant par 
exemple les valeurs de la fonction en ces points. On obtient 
le système 

tip ufti—mpttquptl 
SE OS GRR 0, 


uÿ=Uo(ar), k=1, 2,..., K—1, (18) 
until gti, unti = frs, 


où æ"*1, f"*1 sont donnés. Il est évident que le problème (189) 
donne une approximation du problème différentiel du type 
(177) sur l'intervalle fini x avec les conditions correspon- 
dantes sur ses extrémités. 

Pour démontrer la stabilité du problème aux différences 
(189) il y a licu de l'écrire sous la forme lu = f et de vérifier 
que la solution « est du même ordre de grandeur que f, pour 
1 quelconque. Bien que les équations aux différences (189) 
soient homogènes, il faut ajouter /; dans les deuxièmes 
membres ct montrer que la solution est du même ordre de 
grandeur que f, &, B, Uo (la non-homogénéité apparaîtra 
lors do l'étude de la convergence par suite de l'erreur d'ap- 
proximation). 

Estimons |uf*! |. Supposons que le maximum de cette 
grandeur sur la couche corresponde au point *. Alors, si 
ut > 0, on a 

until — 2unti+untiLO, 
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les valeurs de gauche et de droite étant majorées par la 
moyenne. Si uÿ*! << 0, la dernière expression > 0. En vertu 
de (189) et en ajoutant /;, le signe de cette expression se 
conserve également pour la grandeur 


THÉ — U} _ ja 
T k° 
Par conséquent dans les deux cas 
juntt|<ur + i<lurl+r {fl 
Par hypothèse, dans le premier membre on a la valeur maxi- 


male de |u®*! | sur la (n + Â)-ième couche. En utilisant 
comme d'habitude la désignation 


I] u" [= maxlu |, 


on obtient à partir de la dernière inégalité 
En ESS 8 en 2 EE ESS 22 Emi LE 


Si par contre max [ ut! | correspond à la limite, il est 
égal à at! | ou | Bt! |. Ainsi, 


fuñtt]max(Tartt|, TRI, ut: 


En appliquant cette inégalité pour l'estimation ultérieure 
de uw" à l'aide de u"°1, etc., on obtient 


ur*t]max(llall, BI, 1Z01+(@%+17171) (190) 
avec les notations suivantes 
1/11 = LÉ If" ||, a || — ji) F 


Comme (n + , t LT = const, l'estimation (190) signi- 
fie que le problème aux différences (189) est stable. Remar- 
quons que le résultat obtenu antérieurement à l'aide du 
critère spectral s'est trouvé confirmé. Le schéma aux diffé- 
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rences (189) est stable quelles que soient les relations entre 
les pas t, À. Du point de vue du domaine de dépendance ce 
résultat n'a rien d'étonnant. Pour le calcul de chacun des 
u%'1 il faut résoudre le système d'équations (189) ; ceci étant, 
chacun des «f sera pris formellement en considération. 

Pour d’autres problèmes liés à l'intégration des équations 
d'évolution et des systèmes de telles équations, on peut éga- 
lement trouver des schémas aux différences implicites abso- 
lument stables. La possibilité de choisir les pas t, À con- 
formément à la précision requise permet souvent de réduire 
notablement les calculs, même s'il s'agit de l'augmentation 
de la quantité d'opérations arithmétiques pour chaque 
point de calcul due à la nécessité de résoudre des systèmes 
d'équations. 

Comme nous allons le voir ci-dessous, les particularités 
de ces systèmes permettent lorsque ceux-ci sont linéaires 
d'utiliser des méthodes de résolution efficaces et relative- 
ment simples. Donc lors de l’approximation des équations 
différentielles non linéaires par des schémas aux différences 
implicites ces derniers doivent être linéaires par rapport aux 
grandeurs de la (#7 + 1)-ième couche inconnue. 

Supposons que le problème initial soit quasi linéaire, 
c'est-à-dire linéaire par rapport aux dérivées d'ordre élevé. 
Dans ce cas, en utilisant une approximation implicite pour 
ces dernières, on obtient un problème aux différences dont 
les conditions de stabilité seront déterminées par l'approxi- 
mation explicite des termes et des coefficients d'ordre infé- 
rieur. Généralement ces conditions ne sont pas difficiles 
à remplir. 

Donc si, dans le problème envisagé, le cocfficient de 
conductibilité thermique a est une fonction de U, a = a (U), 
en faisant appel à l'approximation (188) avec a = a (ux), 
on obtient un schéma aux différences implicite, toujours 
linéaire par rapport aux inconnues w%'!. Le problème aux 
différences (188) devient bien entendu non linéaire, tout 
comme l'équation différentielle initiale. En appliquant 
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les méthodes usuelles ($$ 5, 6), on peut faire apparaître que 
la limitation imposée aux pas du réseau, à savoir la condi- 
tion (179), dans ce cas également se trouve levée. 

Il no faut pas penser que l'existence de schémas impli- 
cites rend unutiles les schémas explicites. La simplicité et 
la compacilé de ces derniers sont souvent très précieuses, 
surtout dans le cas des problèmes non linéaires compliqués. 


Problèmes 


Elaborer et étudier différents schémas aux différences 
implicites permettant d'étudier dans l'intervalle 0 Sr <1 
les problèmes suivants: 


LÉ TE +a 0, a>0, 


ot 
U (0, x)}= Ur), U(t, 0) = a (t). 


ou QU Ô oU 
nn Homer «Te 
2. ol je U x Ur U 0x ? 


U (0, x)—Uo(x), Ut, 0)=a(t), U(4, 1)—=B(#). 


17] Ce 
dx 0, 


oU 
d, or Le 
av oÙ 
Dre 
U(0,x)=Uo(x), V(0, x) = Voix), 
U (4, 0) = a (1), V (+, 1) = (+#). 


RÉSOLUTION DES ÉQUATIONS 
$ O. aux Dirrénences 


Lorsque l'on utilise des schémas atix différences implici- 
tes il faut sur chaque couche résoudre un système d'équations 
différentielles. Ce n'est qu'après avoir donné la méthode 
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de leur résolution que l'on peut considérer la description 
de l'algorithme de calcul terminée. 

Le nombre d'équations et d'inconnues est très grand, de 
l'ordre de 1/L. Si l’on ne tient pas compte du caractère spé- 
cifique du système, en le résolvant comme un système du 
type général, un grand nombro d'opérations arithmétiques 
est nécessaire, ce nombre est bien plus grand que dans le 
cas des schémas explicites. On peut montrer que pour les 
systèmes linéaires le nombre d'opérations nécessaires est de 
l'ordre de 1/h°. 

De plus, il ne faut pas oublier qu’un algorithme ne peut 
être réalisé avec une précision absolue, les calculs étant 
toujours menés avec un nombre limité de décimales. Lorsque 
l'ordre du système est élevé, l'accumulation de l'erreur 
d'arrondi peut provoquer une catastrophe. 

Revenons au système linéaire (189) que nous écrivons 
comme suit 


Uo=@, 
— rüpes + (+27) Un — Tunes = U, 
kh=tl,; 25:44, K=1T, 
Ur =. 


Nous avons omis l'indice # + 1 auprès des inconnues et 
introduit la désignation suivante 


(191) 


Comme nous l'avons vu lors de la démonstration de la 
stabilité du problème (189), la solution du système (191) 
doit satisfaire à l'inégalité 


[u|<max(|æ|, |], ques | ur] ). 


Il en découle quo le système homogène correspondant, obtenu 
pour & = f = u} = 0, n’a qu'une solution triviale uy = 0. 
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Par conséquent, la solution du systèmo d'équations (191) exis- 
te et elle est unique. 

Le caractère spécifique du systèmo (191) consiste en 
ce que chaque k-ième équation ne contient que trois incon- 
nues Up 1, Uhs +1 Ceci permet d'exclure successivement w,, 
Us Uo, + . « par la méthode simple suivante. La valeur u, 
est donnée, donc l'équation correspondant à X = 1 ne 
contient en fait que deux inconnues «, et &,, c’est-à-diro 
qu'elle donne la relation existant entre elles. À l'aide de 
cette relation, en utilisant l'équation suivante, on exclut «w; 
et on obtient une relation entre w, et #4, etc. Considérons 
la relation existant ontre uy_1 ot uy: 


Uni = Lattn + An: (192) 


Substituons cette expression dans la Æ-ième équation et ré- 
solvons celle-ci par rapport à w,. On obtient 


Up + u}, + UF 


URL rl 


c'est-à-dire la relation existant entre la paire suivante 
d'inconnues w4 et uy41. Ecrivons-la sous la forme (192) 
posant à cel effet 


: 

Lx = TEL, ? (193) 
u} +rMr 

Me He (194) 


Ces deux formules permettent de passer de ZL4, M, à Ly+, 
M4 pour # quelconque. Comme u, = &, en vertu de (192) 
il y a licu de poser Z, = 0, A1, = et à l’aide des formules 
de récurrence (193), (194) calculer successivement tous les 
Li, M jusqu'à Lx, MX inclus. Puis, comme ux = f, en 
utilisant la formule (192), on trouve successivement tous 
les uz. 

Ainsi, le processus de résolution du système d'équa- 
tions linéaires (191) $o réduit au calcul successif, à l’aide des 
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formules (193), (194), des coefficients Z,, A1, et suivi du 
calcul de «4, à l'aide de la formule (192). C'est pourquoi 
on appelle ce processus méthode du balayage. L'avantage 
essentiel de cette méthode est d'être très économique. Il est 
facile de voir que le nombre d'opérations arithmétiques 
exigé par celte méthode est du même ordre de grandeur 
que Île nombre d’inconnues, soit Æ<1/h, c'est-à-dire est 
minimal. 

Nous allons maintenant vérifier la sensibilité de la 
méthode du balayage aux erreurs d'arrondi, c’est-à-dire 
trouver la relation existant entre la précision des calculs et 
la précision de la solution obtenue. Pour estimer le dévelop- 
pement et l'accumulation de ces erreurs nous allons supposer 
que le calcul approché (avec arrondissement) à l'aide des 
formules (192), (193), (194) peut être interprété comme un 
calcul exact suivant d’autres formules voisines de celles-ci. 

Nous allons commencer par la formule (193) donnant la 
transition de Z,4 à L+,. L'erreur de la valeur cffectivement 
calculée L,+, par rapport à sa valeur exacte apparaît pour 
deux raisons. Primo, la valeur utilisée L, contient l'erreur 
ôLy, provenant des calculs antérieurs. Pour trouver son 
influence sur l'erreur ôL,+,, on prend la dérivée de l’expres- 
sion (193) par rapport à Z4. On oblient : 


ÔLyr = (Lr+) ÔLr. 


Secundo, comme le résultat de chaque opération arithmétique 
est arrondi, même si la valeur utilisée ZL, est exacte, L,4, 
sera entaché de l'erreur apparaissant dans le cycle donné de 
calcul. Désignons cette erreur par Ô4. 

Ainsi, en tous cas pour ÔZL € L,ona 


6x = Li+10Ln + Op. (195) 


La grandeur 6, caractérise la précision des calculs. 

On voit donc que l'évolution de l'erreur d'arrondi se 
trouve entièrement déterminée par les coefficients L,;. En 
particulier, si | L, | > 1, l'augmentation de ôL, est expo- 
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nenticlle, la précision décroissant rapidement. Nous allons 
donner une estimation de la grandeur Z,. 
Comme 7, = 0, en vertu de (193), on a 


r 
Li 2 <1. 


Supposons que 0 << L, << 1. Il est facile de voir, à partir 
de (193), qu'alors on a l'estimation suivanto 


CE La << <1, (196) 


c'est-à-dire que tous les Z, no sont pas supérieurs à r/(1 -- r). 
En utilisant la formule de récurrence (195), on obtient 


2 
ôLu1<lôn1+ (5) Lx 


<|ôw-1|+ (5) 18v-21+ —) ÔLx-21< 


Â r 2N 

Sax | ôn | Fer RE (-—) |ôLo|, 
1+r 

ce qui montre immédiatement que l'erreur ôZ est do l'ordre 


de 6. 
Les formules (194) et (192) pouvent être étudiées d’une 


manière analogue. En dérivant, on obtient 
OM hui = Lan ÔONTs + On, 
OU -, = L,ôur + Ô},; 


où Ô’, Ô” contiennent les erreurs ôu", ôa, ô6fB. En utilisant 
l'inégalité (196), on voit que lors du calcul do A7 et w l'ac- 
cumulation de l'erreur d'arrondi n'a pas lieu non plus. 

Nous avons montré quo lorsqu'on résout le système 
d'équations linéaires (191) par la méthode du balayage, la 
précision du résultat coïncide avec celle des calculs et des 


données initiales. 


101 


Cette affirmation n'est pas vraie pour toutes les méthodes. 
Nous allons voir un exemple. En écrivant chacune des équa- 
tions du système (191) sous la forme suivante 


Un = — Un + SE Un + u?, (197) 


= 
on peut, à partir do la paire de valeurs u,_1, ur, calculer 
Un41. Pour commencer le processus, on a besoin des valeurs 
uo et 1. Pour la première d’entre elles on a &#, — «. La sc- 
conde peut prendre une valeur arbitraire, u, — 0 par exem- 
ple. Les autres u, peuvent être calculées à partir de (197). 
Désignons par uf” la solution obtenue. Il est évident que 
dans la pratique celle no satisfcra pas à la condition limite 
ux —= B. Trouvons la seconde solution u#” de la même 


manière on posant uw!” — 1. Considérons une combinaison 
linéaire de ces solutions 


un = CU + (1 — c) uf? (198) 
qui de toute évidence satisfait à la condition limite à l’ex- 


trémité gauche et aux équations. Nous allons trouver c 
de façon à satisfaire à la condition limite ux = f. On a 


B— u$? 
C—= MEET . (199) 


Alors la formule (198) donne la solution du problème. 

A première vue, la méthode décrite est même plus 
commode que la méthode du balayage. Cependant elle ne 
peut être utilisée ce que nous allons montrer sur un cas 
particulier où l’on peut écrire la solution sous forme expli- 
cite. Posons u# — 0 et cherchons la solution sous la forme 
u, = Const-g". En substituant cette expression dans (197), 
on obtient une équation quadratique en q, soil: 


g—(2++)9+1=0. 


Pour r > 0 quelconque ses racines sont réelles, et co qui 
est important, l’une d'elles est supérieure à l'unité, q <<, 
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{a >> 1. Toute solution du problème est une combinaison 
linéaire de g" et g#, il cest alors facile de montrer que dans 


le cas présent uf” ot uf” mentionnés ci-dessus s'écrivent 


1) aq29} — aq 19% 2) (g2 — 1) g$ + (1 — ag) qù 
D = 
{2 — Qi 2 — Qi 


par conséquent, la solution de (198) est 


Ÿ 


gh À 
ND 2 
Un = U} ny (200) 

Nous allons voir les conséquences des erreurs d'arrondi 
lorsqu'on réalise le processus (197)-(199). Même en idéalisant 
et en supposant que l’on arrondisse seulement la grandeur c 
de (199), tous les autres calculs étant faits avec une préci- 
sion absolue, la solution sera quand même entachée d'une 
erreur égale à 

h_ k 
ÔUr ad Ôc LÉ ben + 9 
I2—q1 

ceci en vertu de (200). Comme g* © 1 pour X grands, une 
petite erreur dans le calcul de c entraîne une grosse erreur 
dans la solution. La grandeur # est de l’ordre de 1/h. Ainsi, 
pour atteindre une certaine précision dans la solution, les 
calculs doivent être menés avec une précision g{/À fois plus 
grande. Par exemple pour qg, = 2 et k = 100 on obtient 
gh — 2100 © 10%, I] est évident que l'on ne peut avoir 
une telle réserve de décimales. C'est pourquoi cette méthode 
correcte du point de vue théorique doit être rejetée. 

Nous avons exposé la méthode du balayage sur l'exemple 
du système (191) correspondant à l'équation de la conducti- 
bilité thermique pour l'intervalle fini (0, X) avec les condi- 
tions aux limites U (4,0) = a (t), U (t, X) = B(t). Pour 
que le lecteur puisse so faire une idée des possibilités de 
généralisation de la méthode exposée à des problèmes plus 
compliqués, nous allons nous arrêter sur certaines questions 
fondamentales. 
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Toute condition aux limites peut être interprétée comme 
l'expression effective de l'influence des processus ayant licu 
dans la partie extéricure, rejetée de l'espace. De co point 
de vuc les équations différentielles décrivent elles-mêmes 
la diffusion de cette influence à l'intérieur du domaine. 
Le problème aux différences est une approximalion du pro- 
blème différentiel. En particulier on peut lui appliquer 
également les considérations sur le rôle des conditions aux 
limites et des équations. La méthode du balayage exprime 
d’une manière explicite ce processus de diffusion de l'in- 
fluence des conditions aux limites à l'intérieur du domaine 
dans le cas où les équations différenticlles sont linéaires. 

Le schéma général de la méthode du balayage pour un 
problème linéaire quelconque est le même que celui qui 
a été envisagé pour les cas de la conductibilité thermique. 
La réalisation aux différences des conditions aux limites 
sur l’une des extrémités do l'intervalle de calcul (par exem- 
ple sur celle de gauche) donne les relations entre les valeurs 
de la fonction discrète (ou du vecteur fonction) aux points 
limite ou voisin de la limite. En utilisant les équations 
aux différences on exclut successivement les inconnues, 
on fait « parcourir » la relation entre les valeurs voisines 
de la fonction discrète à travers tout le domaine de calcul 
jusqu'à l'extrémité droite. Les relations entre les conditions 
aux limites que l'on a ici permettent do trouvor la valeur 
de la fonction au dernier point. En se déplaçant maintenant 
dans le sens inverse, à l'aide des relations obtenues anté- 
rieurement, on trouve successivement toutes les valeurs 
de la fonction discrète. 

Evidemment ce n’est là quo le schéma de la méthode, 
sa forme concrète pouvant changer notablement d'un pro- 
blème à l’autre, suivant les particularités des conditions aux 
limites et des équations. 

Pour la solution des systèmes d'équations aux différences 
on peut également utiliser d’autres méthodes, on particu- 
lier itératives. Revenons au système d'équations (191) et 
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récrivons-le sous la forme suivante: 


Uo—=X, 

r (uns + unes) + ux LD: | 
PARTS Ado k= 1, D'riss K 1 (201) 
Ur — f, 


que nous allons utiliser pour le processus itératif. En substi- 
tuant dans le second membre la v-ième approximation de 


u%), on obtient uWTF dans lo premier membre. Etudions 
la convergence des itérations. losons 


UN = uy + 6), 


où «,, est la solution exacte du système (201). Il est facile 
d'obtenir les formules suivantes: 


65 +100, 
RES EF (6, +61), 1, 2, ..., K—1, 
6%" 0, 


déterminant l'erreur en fonction du numéro de l'itération. 
L'estimation directe de sie donne 
FD v) 
max | 6 |< Fr max | 8 |: 

Par conséquent, Îles itérations convergent, SEE 0, la 
rapidité de convergence . 
O = Lio. TL <1. 

Arrêtons-nous sur le nombre d’itéralions. Supposons que 
l'on prolonge lo processus jusqu'à la coïncidence de deux 


itéraliong successives u$y+1 —— us"), c'est-à-dire que l'on 
obtient la solution du système (201) avec une précision 
maximale. Cependant une telle précision n'est pas indis- 
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pensable car le système (201) lui-même n’est qu'une approxi- 
mation du problème initial différentiel. Lorsque l'on utili- 
se la méthode du balayage on obtient cette précision gra- 
tuitement, ici au contraire pour l'obtenir il faut faire appel 
à des itérations successives, par un travail supplémentaire 
pouvant par conséquent être superflu. 

D'un autre côté, la convergence de la méthode numérique 
(découlant de l’approximation et de la stabilité) a été démon- 
trée en supposant que le problème aux différences était exact. 
II est facile de montrer que si l’on tient compte des erreurs 
d'arrondi, cela ne change rien (car la stabilité les amortit). 
L'influence des erreurs apparaissant par suite des itérations 
insuffisantes doit être étudiée spécialement. 

En exagérant, supposons que l'on so borne à une seule 
itération en prenant pour l’approximation initiale uf° = ux. 
Ceci signifie qu'en fait pour le calcul do u?‘!, on utilise 
la formule suivante [voir (201)]: 


r(uy_ Huy) +ur 
14 2r | 
En estimant [ut] on obtient 


(202) 


n+1 — 
U, — 


maxfuñ+i|<max]un], 
k h 


c'est-à-dire le schéma aux différences (202) est stable (pour 
r quelconque, bien que le schéma soit explicitel). 

Revenons à l’approximation. En substituant dans (202) 
les développements suivants 


Ur" =U+U+0(r), 


k? 3 
Ua = UÙ EU ht + Ur Æ Usa + O (R'), 
on obtient 
aU xx + 0 (h?) 
1+2r ° 
Par conséquent, la relation aux différences (202) ne donno 
une approximation de l'équation de la conductibilité thermi- 


Ui+0 (x) — 
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que que pour r —+ 0. Sir par contre est fini, en utilisant (202) 
on aura une autre équation. On peut montrer (voir pro- 
blème 4) qu'il en sera de même après la v-ième itération, 
c'est-à-dire que l'erreur d'approximation est une fonction 
de ret de vet ne tend vers zéro que pour r — 0 ou v —+ oo. 
Ainsi, un nombre insuffisant d'itérations pout donner lieu 
à l'absence d'approximation. 

Un autre cas est possible. Supposons, par exemple, que 
pour la solution du problème (191) au lieu de (201) on ait un 
autro processus itératif du type 


u(s+1) — (p (un ufs), uv) un), 


1” h+1? 


pour lequel la condition d’approximalion se trouve vérifiée. 
Il est évident que dans le cas présent on risque de ne pas 
vérifier la condition de stabilité. Pour cstimer le nombre 
nécessaire d'itérations nous allons raisonner do la manière 
suivante. Lors de la première itération, pour trouver uf? 
on tient compte en plus de u% seulement do uÿ.1 = uf1 
et de uk, — u}®%,. Lors do la seconde, en utilisant w}?: 
et uw, on tient compte par là même do uf. = uf”, cet 
Ufta = Uf», Cte. Chaque itération étend le domaine de 
dépendance d'un point de calcul, ceci d’un côté et de l’autre. 
C'est pourquoi v itérations sont équivalentes à un certain 
schéma explicite utilisant 2v + 1 valeurs de la #-ième 
couche. Pour la stabilité de ce dernier, il faut vérifier une 
certaine inégalité du type 


T 
—— < const 
(vh)2 < . 
c'est-à-dire r & const -v*. Ainsi le nombre nécessaire d'itéra- 
tions est égal à 
v mr. 


On voit que lorsque r n’est pas trop grand, les méthodes 
itératives peuvent donner des résultats satisfaisants et même 
faire concurrence à la méthode du balayage car elles per- 


107 


mettent de diviser le processus de calcul en plusicurs bran- 
ches parallèles ct n'ont pas besoin d'une grande mémoire 
fixant des massifs importants de cocfficients Z, AJ indispen- 
sables pour la méthode du balayage. Il est parfois commode 
de combiner les deux méthodes en utilisant le balayage sur 
les domaines isolés et de les « ajuster » par itérations. 

Nous allons faire une dernière remarque. Nous n'avons 
étudié que les équations aux différences linéaires. Les métho- 
des itératives ont un champ d'application plus large, par 
contre la méthode du balayage utilise en fait la linéarité 
des équations. Pour la résolution des équations non linéai- 
res seules les méthodes itératives conviennent. Cependant 
les itérations peuvent être réalisées différemment. Ci-dessus 
nous avons envisagé un algorithme explicite du type 
u(W+D = p(uM). Pour r grand il converge lentement car 
il faut tenir compte du domaine de dépendance de la solu- 
tion. Ce dernier est donné par les propriétés du problème, 
celles-ci apparaissant déjà dans son modèle linéaire. Ayant 
à notro disposition un instrument aussi efficaco que la 
méthode du balayage, on peut lors de l'élaboration du pro- 
cessus itératif utiliser des algorithmes implicites linéaires 
et ramener ainsi le problème à la résolution d'un système 
linéaire sur chaque itération. Ceci permet d'obtenir un 
processus à convergence rapide. 

Supposons, par exemple, que dans le problème envisagé 
ci-dessus le coefficient a dépende de la fonction cherchée 
a — a (U). Il est évident que lors de l'approximation par un 
problème aux différences on peut prendre a (u"), mais 
supposons que pour des raisons quelconques ceci ne nous 
convienne pas et que l'on utilise a (u"+1). Le système d'équa- 


tions aux différences (191) se trouve alors êtro non linéaire 


r = ru). Il est tout naturel de poser r == r (uÿ”) ot on 


résolvant le système (191) par la méthode du balayage 


d'obtenir ufÿ*), Si pour l’approximation initiale on prend 
uf” = uf, la rapidité de convergonce du processus itératif 
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UV) + uf# sera déterminée par la grandeur do la variation 
effective de la solution d'une coucho à l'autre. 


Problèmes 


1. Estimer l'ordre de grandeur du nombre d'opérations 
arithmétiques nécessaires à la résolution du système de # 
équations linéaires du type général par la méthode d’exclu- 
sion. 

2. Si la méthode do résolution d’un système d'équations 
obtenu à partir d’un schéma implicite est donnée, l'opéra- 
teur de transition d'une couche à l'autre se trouve ainsi 
déterminé, on peut alors étudior le schéma en tant qu’expli- 
cite et l'écrire comme suit 

unti= Nu" +rf"). 
En utilisant la stabilité du problème, c'est-à-dire la condi- 
tion limitant les puissances de l'opérateur /?, [| R" || < const, 
étudier le comportement des erreurs apparaissant dans u"*1 
par suite de l'arrondissement du second membre lors de 
l'augmentation de n, nt <t. 

3. Etudier l'équation aux différences 

unti= qi, 


n+i n n+1 n+!i 
Uy, — y se Unie VR 1/9 —( 
T h FA. 9 

TU —v} ut} —_unti k=1, 2, #7 K, 
k-—1/2 k—1/9 k R-1 =0 
— V9 


T h 
N+1  __ an+1! 
UK +172 =D 
OÙ Uyy17, Signifie que cette valeur correspond au point 


Tati = (% + 1/2) h. Appliquer la méthodo du balayage 
en utilisant une relation du type 


un = Livni + M, où ui = Lux + As. 
Vérifier si l'algorithme de calcul est correct par rapport 
aux erreurs d'arrondi. 
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Trouver le schéma analogue pour la résolution du pro- 
blème non linéaire suivant 


oU OP (1 
+0,  U(O, 2)=Uo(x), 
ss =0,  V(0,x)=V 
D + = Û ( ? x) = o(t) 
sur l'intervalle 0 < x < X sur les bords duquel on a les 


conditions aux limites suivantes 
U ({, 0) = a (1), V (4, À) = B (4). 


4. Supposons que, lors de la réalisation du processus 
étudié ci-dessus (201) 
r (u (9) 4 + uf 1 +uÿ) 
DT par Fons 


on se borne à W'itérations, c'est-à-dire "+1 = uN) nf = nf 
(pour plus de simplicité omettons les conditions aux 
limites). 

Vérifier que le schéma aux différences obtenu est stable 
pour r ct Ÿ quelconques. 

Montrer que l'erreur d'approximation est proportionnelle à 


or N 
1 —2r | 
A cet effet poser 
9 = a + hb0 5 COR. 
où de toute évidence 
av=ur, DOS Qu), CO (uxx)n, 


et trouver a(N) à l’aide de la formule du processus itératif. 
Comme 


v+1 
ufr) 


until = QU = un + T (u)r +... 


connaissant a) il est facile de trouver l'erreur d'approxi- 
mation. 


Chapitre III 


CALCUL DES SOLUTIONS 
$ 10. pisconrinurs 

Dans notre étude antéricure nous avons supposé que la 
solution exacte du problème inilial est une fonction régulière. 
Lors de l'étude de l'approximation et de l'élaboration des 
modèles linéaires nous avons souvent utilisé cette hypothèse. 
Pour les problèmes différentiels c’est tout naturel, la solu- 
tion devant être telle qu'au moins les dérivées figurant dans 
l'équation existent. Si la solution contient des singularités 
la rendant non régulière, il y a lieu de procéder à une étude 
spéciale et en cas de nécessité de modifier la méthode. On ne 
peut rien dire de plus précis, chaque singularité ayant ses 
particularités. 

Ici nous nous bornerons à l'étude de cette question dans 
lo, cas important pour la pratique où la solution est une 
fonction lisse par morceaux avec discontinuités. Considé- 
rons l'exemple 


HU =0, U(0,2)=Uo(x) (203) 


Dans le cas où la solution U ({, x) est une fonction régulière, 
on est ramené au problème que nous avons déjà étudié en 
détail. Cependant dans certains cas les conditions de régula- 
rité et d'existence de la solution peuvent s'exclure mutuelle- 
ment. Ainsi la solution du problème (203) satisfait évidem- 


ment au système 


dx du 
ee 


c'est-à-dire qu'elle est constante le long des droites 
z = 2 + Uo (to) t. 
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Fig. 12 


Ces lignes sont appelées caractéristiques. Le long de ces 
caractéristiques les équations du problème dégénèrent en 
relations entre les différentielles de la fonction (ou des fonc- 
tions). Dans notre problème la pente de chacune des caracté- 
ristiques est donnée par la valeur de la fonction U, (x) au 
point æx,, point d'intersection d'une caractéristique considé- 
rée avec la ligne des données iniliales { — 0 (fig. 12). Il est 
facile de voir que dans le cas où U, (x) est, ne serait-ce que 
sur un petit segment do l’axe x, une fonclion décroissante, 
les caractéristiques issues des points de co segment sc ren- 
contrent. Chacune d'elles donnant indépendamment des autres 
une valeur de la fonction, la solulion au point d'intersection 
no sera pas univoque. Pour parer à ces difficultés, pour les 
problèmes décrivant des processus physiques réels (hydro- 
dynamique), on admet l'existence do solutions discontinues. 

En particulier la non-unicité mentionnéo indiquera l’ap- 
parition d'une discontinuité. Mais si l’on suppose que la 
solution contient une discontinuité, au point de disconti- 
nuité les dérivées ne sont pas déterminées et par conséquent 
les équations différentielles perdent leur sens. Nous devons 
alors remplacer ces dernières par les relations finies reliant 
les valeurs des fonctions de part et d'autre de la disconti- 
nuité. En prolongeant l'analogie existant avec les problèmes 
hydrodynamiques, il y a lieu de supposer que ces relations 
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doivent exprimer, pour les solutions discontinues, les mêmes 
relations physiques que des équations différentielles pour 
les solutions lisses. Du point de vue mathématique formel 
ces solutions s'obticnnent de la manière suivante. 

Soit U(x,t) une fonction lisse satisfaisant dans un 
certain domaine du plan x, £ à l'équation (203). Intégrons 
(203) sur une partie quelconque $S de ce domaine. Il est 
facile d'obtenir 


J] (+0 GE) dr di = J] (+ (5) ) = 


= ( | dt dx + | | _. (TS) dat Udr—© dt, 
8 8 r 


où la dernière intégrale curviligne est prise sur le contour F, 
frontière de S. Par conséquent si U (x, ?) est la solution de 
(203), on a 


4 U dx — T dt —0 (204) 
T 


pour un contour quelconque l. Lorsque l'on utilise (204) 
au lieu de (203), la fonction U (x, f) n’a pas besoin d'être 
ni lisse, ni continue. C’est pourquoi il y a lieu de l'utiliser 
pour obtenir des relations sur la discontinuité. 

Nous allons raisonner comme suit. La discontinuité dans 
la solution une fois apparue se déplace ensuite décrivant 
dans le plan x, { une certaine courbe z=X ({). Considérons un 
petit élément de cette courbe et construisons sur cet élément 
comme diagonale le contour rectangulaire l (fig. 13). Comme 
ce rectangle est petit, sur chacune de ces deux moitiés on peut 
supposer que la fonction U est constante et égale à U- à gauche 
et U* à droite. Calculons l'intégrale (204) pour ce contour, 
on obtient 


U*Ar—(+)" AU ar+() At=0, 
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Fig. 13 


où Az, At sont les côtés du rectangle. Lorsque l’on contracte 
le rectangle en un point, le rapport Ax/At tend vers X”(#), 
à la limite la dernière égalité donne 


2 2 \— 
wux-() (Æ). eo 

Comme nous l'avons vu ci-dessus, l'intersection des carac- 
téristiques et l'apparition de la discontinuité n’ont lieu que 
si U-> U*. On peut montrer que (pour notre problème) 
seules des discontinuités de ce type peuvent exister. En sim- 
plifiant (205) par U*— U- et en y ajoutant l'inégalité 
mentionnée, on obtient 


X' = U-> U*. (206) 


Cette relation entre U-, U* et X” remplace sur la ligne de 
discontinuité x — X (t) l'équation différentielle (203). 
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Généralisons le problème (203) ct supposons qu'il y ait 
des discontinuités sur certaines lignes x — X (t) satisfaisant 
aux relations (206). 

Nous allons passer à l'élaboration de la méthode numéri- 
que de résolution des problèmes (203), (206). Le plus simple 
est d'utiliser les méthodes étudiées ci-dessus en les modifiant 
soulement au voisinage immédiat de la ligno de discon- 
tinuité. Ce n'est pas difficile, tout au moins pour notre 
problème simple, cependant un certain nombre d'incon- 
vénients apparaissent. Los formules spéciales pour le calcul 
des grandeurs sur la ligne de discontinuité et dans les points 
voisins doivent tenir compte de toutes les dispositions 
possibles de cette ligne par rapport aux points du réseau 
de base. Do plus il faut envisager l'apparition éventuelle 
d'une discontinuité et par conséquent vérifier la solution 
obtenue d'une façon convenable (dans tous les points de 
calcul). Il est évident que cela augmente notablement le 
volume de l'algorithme de calcul qui perd alors sa simplicité 
ot sa compacité, vu le nombre restreint de points de calcul. 
C'est pourquoi souvent on préfère une autre méthode pour 
obtenir les formules de calcul que nous allons décrire ci- 
dessous. 

Pour le problème aux différences la notion de discon- 
tinuité de la solution du point de vue formel n’a pas de sens, 
la fonction discrète étant déterminée sur un ensemble discret 
de points, le réseau de calcul. D'un autre côté, comme nous 
l'avons déjà vu, la condition sur la discontinuité (206) dé- 
coule de l'équation différentielle (203) et par conséquent y 
est « contenue ». La théorie des équations différentielles nous 
apprend que l’on peut obtenir la solution discontinue comme 
la limite de la solution lisse de l'équation perturbéo lorsque 
le paramètre perturbateur tend vers zéro. Nous pouvons uti- 
liser ce fait car en utilisant telle ou telle méthode numéri- 
que on remplace toujours le problème initial par un autre 
problème discontinu perturbé. Nous allons réaliser l'ap- 
proximation en deux étapes. Tout d’abord en introduisant 
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une perturbation, on remplace le problème initial par un 
problème intermédiaire, puis on passe au problème aux 
différences. 

Au lieu de a considérons l'équation 


Davis (Lo eo 


où e est un paramètre petit. I] ost évident que si l'on so 
borne à des fonctions lisses, alors, & étant petit, les solutions 
des problèmes (207) et (203) pour des conditions initiales 
identiques ou voisines seront également voisines. Pour se 
faire une idée do la différence apparaissant dans le cas des 
solutions discontinues (pour (203)), considérons le cas parti- 
culier suivant. 

Supposons que la solution du problème (203), (206) soit 
une fonction en escalier 


UT pour z—ot<o0, 
9 ea pour æ—>0, (21e) 
où 
= EU y-> ur, |(209) 


2 


et U-, U* sont des constantes. Il est évident que la fonction 


(208) satisfait à (203), (206). 
Cherchons la solution correspondante de l'équation (207) 


sous la forme 


U,(x,t) = { (x — wt). (210) 
Il est naturel d'admettre que 
f(x) U+t pour zx—>-—+oo, (211) 


car loin de la discontinuité les solutions U, U, sont des 
fonctions lisses et par conséquent voisines. En substituant 
(210) dans (207), on obtient une équation différentielle 
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ordinaire pour / 
— of" + ff +e (7) = 0. (212) 


Il est facile de voir que 


; x —(d£ 
= @+const.sin — 
{ = © + V2 


est la solution de (212). Comme 
f = const 


salisfait également à l'équation (212), la solution qui nous 
intéresse est de la forme (fig. 14): 


— (DE TT 
U” jour _—— 
l e 1/2 2 ? 
Ut+U- Ut—U- . xz—ot x — wt ELA 
——— + —— sin — pour | —— — 
U, — 2 + 2 e V2 P Gé < 2 ? 
U* | pour 75 >+. 


Cette solution est une fonction lisse, au lieu de la discon- 
tinuité U-, U* elle a une zone de transition continue de 


U- à U* dont la largeur est ex V 2. Si e est suffisamment 
petit, ce domaine est étroit et U, cest voisin de la solution 
discontinue (208) du problème initial (203), (206). 

Ceci permet de remplacer les problèmes (203), (206) par 
le problème (207). Il est très important que la solution de 
ce dernier est une’ fonction lisse. C'est pourquoi lors de 
l'élaboration de la méthode numérique on peut utiliser 
toutes les méthodes étudiées ci-dessus permettant de trouver 
et d'étudier les équations aux différences. L'introduction 
d’un terme perturbateur (appelé viscosité artificielle) compli- 
que les calculs mais on revanche donne la possibilité d'ef- 
fectuer le calcul à l’aide de formules bion connues sans faire 
appel à l'étude des particularités et aux épreuves sur l’appa- 
rition des discontinuités. 
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AE x=-wt 
Va Vz 
Fig. 14 


La méthode envisagée faisant ressortir la viscosité arti- 
ficielle (207) (qui n'est pas la seule possible) est commode 
parce que la largeur de la zone « floue » de la discontinuité 
est de l’ordre de e et ne dépend pas de la grandeur de la 
discontinuité U- — U+. Les variations de la fonction dis- 
crète sur des intervalles supérieurs à æe ont alors un sens 
réel pour une solution quelconque. C’est donc sur ces considé- 
rations que doit être basé le choix de la grandeur e. 

Pour ce qui est du choix de k (et + pour des schémas 
implicites) remarquons que l'introduction de la viscosité ne 
donne l'effet désiré que dans le cas où la zone de disconti- 
nuité contient au moins plusieurs points de calcul. Dans 
le cas contraire il n'y a pas lieu d'attendre une approxima- 
tion convenable, l'équation aux différences pouvant réagir 
non correctement à la viscosité. 

Ainsi on a une certaine relation entre les paramètres e, 
%, k, c'est-à-dire e — e (x, k). D'un autre côté tout problème 
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aux différences a une certaine viscosité (dite d'approxima- 
tion) car le schéma aux différences est équivalent aux équa- 
tions initiales majorées de l'erreur d'approximation, par 
exemple © (rt, k). Cette dernière est en réalité une certaine 
expression différentielle à coefficients petits, c'est-à-dire 
elle peut être interprétéo comme une viscosité artificielle. 
C'est pourquoi, si l'on veut, la méthode exposée ci-dessus 
d'élaboration du problème aux différences se ramène à la 
méthode générale avec une erreur d'approximation d'un 
type spécial. 

Nous avons étudié un seul exemple (203) particulière- 
ment simple, néanmoins nous avons exposé tout ce qui est 
important sur les méthodes de solution de ce type de pro- 
blèmes lorsqu'il y a des discontinuités, à une exception 
près. 

Bien que ça puisse paraître étrange, ce sont les équations 
linéaires ct en général les problèmes où la discontinuité 
s'étend le long des caractéristiques (c'est-à-dire la ligne 
de discontinuité x = X (f) est une caractéristique) qui font 
exception. Dans ce cas les discontinuités ne peuvent être 
étalées d'une manière stable à l'aide d'une viscosité arti- 
ficielle. C'est pourquoi en cas de nécessité, pour le calcul 
correct d’une telle discontinuité il y a lieu d'utiliser des 
formules spéciales. 

Pour le calcul d'une singularité quelconque deux possibi- 
lités s'offrent. Ou bien on donne une description détaillée 
de la discontinuité ou bien on l'étale. Dans ce paragraphe 
on a étudié un exemple du second type. 


Problèmes 


1. Pour les problèmes (203), (206) trouver les formules 
aux différences pour le calcul des grandeurs sur la disconti- 
nuité et aux points voisins. Une maille de calcul typique est 
donnée sur la fig. 15. Les points de calcul disposés sur la 
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nt 


Fig, 15 


ligne de discontinuité sont désignés par des croix. On y 


calcule deux valeurs, de gaucho w- et de droite w 
2. Pour le calcul E” solutions discontinues de l'équation 


nn 
T + SE =0 
avec 7 condilions suivantes sur la ligne de discontinuité 
z = X (t) 
(U*—U") À'=F(U*)—F (07, 
Fo(U)> Fu (U”), 
où F est une fonction donnée U, appliquer les différentes 
méthodes possibles d'introduction de la viscosité artifi- 
cielle. En plus de 
aU , 0F CU a f aU \2. 
on ge de ge \i0r ) — 


ot 
envisager l'équation Re du type 
OU , 0F(U) _. &U 
FT 'RRRAT EE 
Dans les deux cas pour la fonction U, (x — owt) étudier les 
conséquences de l'introduction de la viscosité pour la solu- 
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tion discontinue. Etudier séparément les mêmes questions 
pour F (U) linéaire, F = aU. 

3. Trouver le schéma aux différences pour l'équation 
(207). Etudier l'approximation et la stabilité. 


$ 11. PROBLÈMES MULTIDIMENSIONN ELS 


Lorsque l'on passe des équations différonticlles ordinai- 
res aux équations aux dérivées parlielles, c'est-à-dire lorsque 
l'on passe d'une variable indépendante à deux, comme nous 
l'avons déjà vu, non seulement des difficultés quantitatives 
apparaissent, mais également de nouveaux problèmes im- 
portants. Les problèmes essentiels ont été étudiés dans les 
paragraphes précédents. Le passage aux problèmes à trois 
et plus variables indépendantes n'est pas non plus trivial. 
Cependant presque toutes les méthodes envisagées d'élabora- 
tion et d'étude des équations aux différences peuvent être 
simplement et naturellement généralisées à ce cas. Les 
problèmes à deux variables indépendantes #, x sont à cet effet 
un bon modèle de problème multidimensionnel. 

Nous allons brièvement nous arrêter sur les questions 
essentielles apparaissant dans les problèmes où les variables 
indépendantes sont le temps { et deux coordonnées spatiales 
z, y. Ces problèmes sont appelés bidimensionnels. Pour la 
démonstration nous allons prendre l'équation de la conduc- 
tibilité thermique, c'est-à-dire le problème 


Ut SU , 2U 
GE GT to ? U(0, x, y)=U(x, y). (213) 


Dans le cas bidimensionnel le réseau de calcul simple 
se composera de points de coordonnées 
Ent, Zn=khx, Ym=mMhy, 
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n,t 


Fig. 16 


il est, par conséquent, déterminé par trois paramètres 
T,h;, k, qui sont les pas du réseau. Désignons par 4, » 
les valeurs de la fonction discrète correspondant à ces points. 

Toutes les méthodes énumérées au $ 7 d'élaboration des 
formules de calcul peuvent directement être généralisées au 
cas multidimensionnel. 

En particulier, pour le problème (213) lorsque l’on uti- 
lise la maille de calcul représentée fig. 16, n'importe laquelle 
de ces méthodes donne l'équation aux différences 


n+i n 
Um Uk, m 
T 
n n n 
=. Up, 4.1 | m2, nt Up _1, m | uy, m+i 2uÿ, m + up, m-i{ (214) 


h2 h} 
0 
Uk, m — Us (za, Um)» 
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qui est une généralisation directe du cas unidimension- 
nel (178). 

Le schéma de principe de l'étude de la convergence rame- 
nant ce problème à l’approximation et à la stabilité a été 
exposé dans le $ 5 sous une forme très générale, de telle 
sorte que maintenant le problème bidimensionnel s’en trouvo 
être un cas particulier. [1 suffit de remplacer les deux para- 
môtres Tt, L et les deux arguments de la fonction discrète 
n, k par les trois paramètres T, ,, L, et les trois arguments 
n, k, m. 

Pour la vérification de l’approximation du problème aux 
différences et du problème différentiel on procède d'une 
manière analogue. Ainsi pour les problèmes (213), (214) 
supposant la solution exacte U (t, x, y) lisse, on peut écrire 


Ut = Um +7 (5) +0 (x), 


k, 


n 


Us, m = UF, "ét ( ôx }, m 


+ 
4e (RL) ae (HE) oi 


UF, m£1 = U}, m + hy($ dy a + 


| o2U 1 GU \n 
+ (re (Si + 0 Gb. 


En substituant ces expressions dans (214) on trouve facile- 
mont qu'elle est satisfaite à © (rt, k£, k}) près, c'est-à-dire 
qu'il y a approximation. 

En ce qui concerne l'étude de la stabilité, on fait, nous 
n'avons étudié (au $ 6) qu’une méthode générale, à savoir 
le critère spectral de stabilité des équations aux différences 
linéaires de structure stratifiée du type 


unti — = Ru" + tf". 
123 


Si par fonction sur la couche uw" on entend maintenant 
UR m, C'est-à-dire la fonction discrète de deux indices k, m, 
tous les raisonnements faits au début du $ 6 restent vrais 
et la stabilité sera une borne pour les normes des puissances 
de l'opérateur À. Pour les opérateurs À correspondant à la 
formule 


(Aux, m — À Xp, qUhtp, mtq (215) 


qui est une généralisation de (123), on peut. trouver l'estima- 
tion de ces normes à l’aide du rayon spectral des opérateurs. 
Ainsi il est facile de voir que les fonctions discrètes 


Un, m = Uo, 6 RP +M) (216) 


pour ®, Ÿ quelconques sont les fonctions propres de l'opéra- 
teur /? (215) et 
A = D ap, ae (PP+a), 
P, q 


les valeurs propres correspondantes. L'inégalité | À (p,%) | < 
< 1 donne la condition nécessaire de stabilité. Ainsi dans 
le cas bidimensionnel, le critère spectral de stabilité con- 
serve son efficacité. 

Appliquée au problème (214) cette procédure donne la 
condition de stabilité suivante 


1 
RTS 


L'approximation utilisant les équations aux différences 
implicites donne alors des schémas stables quelles que soient 
les relations entre les pas des réseaux. Pour le problème (213) 
on aura de toute évidence le schéma suivant (fig. 17) 


n+i n +1 n+i ne 
Uh, m— Un, m _ THE m— 2Uk, MT ur .n se 
L h2 
u" — ut +unt 
+ h, ii h,m L (217) 


hÿ 
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noel 
Up _),m 


m.y 


k,x 
Fig. 17 


Sa stabilité pout être vérifiée directement par l'estimation 
de u"*!, comme dans le 8 8. Si pour u" on prend uno fonc- 
tion du type (216), il est facile de voir que u"*! — Au", avec 


A 
L 

4 sin? TT sin? + 
ot pour Tv, L,, k, quelconques on à |A | < 1. 

Ainsi les questions essentielles liées à l'élaboration et 
à l'étude des équations aux différences ne changent pas en 
principe lorsque l'on augmente le nombre de dimensions. 
Néanmoins l'augmentation du volume du problème et la 
complication des algorithmes de calcul peuvent donner lieu 
à de nouveaux problèmes. 

Nous allons nous arrêter à cet effet sur les méthodes do 
résolution des systèmes d'équations aux différences apparais- 


125 


Fig. 18 
sant lorsque l'on utilise les schémas implicites. Considérons 
les équations (217) dans le rectangle (fig. 18) 

1<k<K, 1<m<M, (218) 


en se donnant les valeurs de #"*! sur ses bords. Pour simpli- 
fier on pose À, — k, = h, on introduit la désignation r — 
= t/h?, on omet l'indice # + 1 auprès des inconnues wñ*}, 
et on peut alors écrire ces équations comme suil 


 'U, m4 —7Uh-1, m + (1 + 4r) Um lU4t, m — 
— Tüy,m-1 = Uñ, me (219) 
Les valeurs aux limites de uy,, , sont données comme suit 


Uh, Af+1 — Ô, 


Uo,m— mr Uk+i, m = Pms (220) 
Un, 0 — Yk: 
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Les indices X, m prennent l’ensemble de valeurs (218). Ainsi 
le nombre d'équations et le nombre d'inconnues sont égaux 
à XM, à chaque point intérieur du rectangle correspond 
une équation propre (219). 

Dans le cas unidimensionnel, pour la solution d'un sys- 
tème analogue d'équations on a pu utiliser la méthode efficace 
d'exclusion, à savoir la méthode du balayage ($ 9), car nous 
avions un système très simple, où chaque équation ne con- 
tenant que trois inconnues de numéros successifs, corres- 
pondant à lour disposition naturelle. Ici il n'en est rien, 
néanmoins on peut élaborer une méthode de résolution du 
système d'équations (219) généralisant la méthode du 
balayage qui nous a si bien servi dans le cas unidimen- 
sionnel. 

Nous allons considérer l'ensemble des valeurs uy, :, 
Uh, 91 + +, Up, M pour k donné, comme les composantes 
du vecteur w, de dimension 7 (fig. 18). Choisissons parmi 
toutes les équations du système (219) celles correspondant 
à cette valeur de #. Il est évident qu'elles relicront les com- 
posantes de trois vecteurs seulement u} _;, un, Un+,. Ecrivons 
ces A1 équations sous la forme d'une seule équation vecto- 
rielle 


— Aux + Burn — Cuy,s = dx, (221) 


où À, D, C sont des matrices carrées, et d,, un vecteur 
d'ordre À. 
Il est évident que À = C = cZ (I élant la matrice unité), 
et 
(A+ér —r 0 0) 
—r À+4r —7r., 
(0 —r 1+4r:. 


Le | 1 4r 7 
(0 —r 14+-4r) 
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(us, + rya 
Uh, 2 
U, 3 


dr 


Un. M-1 
ur. M +70 


Comme le système (219) s’est réduit à ÆÀ équations (221) 
reliant des vecteurs uy_1, Un, Un+1 par trois, on peut utiliser 
la méthode du balayage, tenant compte évidemment du 
caractère vectoriel des équations (221). 

Supposons qu'entre up et #y on ait la relation 


Un = Llg + MR, (222) 


où Ly est une matrice carrée, ot M, un vecteur du même 
ordre que w,. En substituant (222) dans (221), on exclut 
ur et l’on obtient donc une relation entre la paire suivante 
de vecteurs 


(B + AL) Up — Cu = A, + di. 


En résolvant cette dernière par rapport à w,, c'est-à-dire en 
multipliant à gauche par la matrice inverse de B — AL,, 
on obtient 


un =(B — AL)! (Cuyss + AM + da). 


Pour écrire cette dernière relation sous la forme (222), 

on pose 

Lau =(B—AL,Y!C 

kH = ( x) . | (223) 

Nas == (B — À Ly) (A My + dh). 
La méthode de résolution est maintenant claire. La con- 

dition aux limites pour À = 0 détermine ZL; = 0, M, = «&. 

A l'aide des formules (223) on trouve successivement tous 

les L,, Af,. Comme on connaît ux+1, Ua = BP, à partir 

de Ly, My en utilisant la formule (222) on obtient toutes 
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les solutions &, du système. Les coefficients L, étant des 
matrices, la méthode exposée est appelée méthode du balayage 
matriciel. 

Apparemment elle est analogue à la méthode unidimen- 
sionnelle du balayage. Cependant à l’opposé de celle-ci elle 
est très rarement utilisée, ce qui est dû à ce qu'elle conduit 
à des calculs compliqués. Lors de chaque cycle du processus 
de calcul il faut prendre l'inverse d'une matrice d'ordre 
élevé —<1/h (et fixer dans la mémoire <1/h de ces matrices). 
C'est pourquoi souvent l’utilisation des méthodes itéra- 
tives de résolution des systèmes (et même l'utilisation des 
schémas explicites) se trouve être plus efficace. 

Lors de l'élaboration d'un algorithme plus économique 
pour le type de problèmes envisagés on procède différem- 
ment. Considérons le schéma aux différences (fig. 19) 
Un — Un, m URSS, ma RM UR, m 
a 12 + 

lx 


+ uR m+i —2uy, m+ux. m1 (224) 

h£ | 
qui est intermédiaire entre (214), (217) donnant également 
une approximation de l'équation différentielle (213). En 
substituant dans (224) les fonctions discrètes du type (216), 
on trouve u"tl — Au", où 


ne sin? RUE 

h 2 
À = Te . 

AT oo P 

1+ F3 sin p 


et par conséquent le schéma aux différences (224) ne peut 
être stable que si 


1/2 9—01279 129 


n,t 


Fig. 19 


ce qui n’est pas étonnant car ce schéma est explicite par 
rapport à l'une des variables, à savoir y. On utilise de tels 
schémas si le problème impose que les directions x et y ne 
soient pas équivalentes. Par exemple si la solution dépend 
faiblement de y, on peut prendre À, bien plus grand que À, 
la condition (225) n'est alors pas difficile à remplir. 

Remarquons le fait que d'un côté le schéma (224) n’im- 
pose aucune restriction à la relation existant entre + et k. et 
d'un autre côté, l'ensemble des équations (224) pour chaque 
m donné forme un système pouvant être résolu par la méthode 
simple du balayage. 

On peut dire ainsi que le problème de l'élaboration d’un 
algorithme efficace est à moitié résolu. Il ne reste plus que 
la condition imposée à v et k, (225). Le nombre d'opéra- 
tions arithmétiques nécessaires pour obtenir la solution se 
trouve être proportionnel au nombre de points de calcul. 

Changeons de rôle les directions x et y, c'est-à-dire en- 
visageons un schéma explicite par rapport à x et implicite 
par rapport à y. On obtient un schéma permettant de résoudre 
la seconde moitié du problème. Les schémas s'excluent 
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mutuellement, cependant nous allons essayer de les utiliser 
successivement, l’un pour les pas pairs et l’autre pour ceux 
impairs par rapport à {. Comme le cycle élémentaire sera 
dans ce cas une paire de pas, il est plus commode d'appeler 
un pas dans le temps tout le cycle, d'utiliser chaque pas pour 
un déplacement de v/2. La solution obtenue sur la première 
moitié du pas est alors une certaine solution intermédiaire u. 
Ceci donne le schéma aux différences suivant 


Fa n ph Condé 
UR, m—Up, m __ Uh:4, m— Jun, m+Uk-1, m + 


T/2 h? 
un — Au} LL +ur 
+ hk,m+1 Am k, mi , (226) 
h;, 
Rs _ us 2e 
Um Uh mi _Ukss, m— 2UR, m7 Uh-1, m + 
T/2 ou h2 
unti —9Junti un+i 
FA hk, m+1 km TU, m—i x (227) 


hÿ 
Comme précédemment, c'est une approximation de l'équa- 
tion (213). Nous allons étudier la stabilité de ce schéma. 

En utilisant pour «" une fonction du type (216), on ob- 
tient à partir de (226) u = Au", où 


2T . , Ÿ 

PR SU LES 

. 1 HS sin D 
: ji sin? L 
hà 2 


et à partir de (227) u*t = A, où 


2T () 

rer 21 

7 _° PE sin o 
2 sim Ÿ 
h3 2 
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Seul le produit ÀÀ = À nous intéresse car c’est lui qui cor- 
respond à un pas entier dans le temps. Il est facile de voir 
que À est le produit de deux facteurs 


2T . , 2T . » 
| RE sin p { 12 sin D 


D 
2T . » ® 2T . , Ÿ ! 
1+ DE sin 1+ DE] sin D 


(228) 


2 


chacun de ces facteurs n'étant pas supérieur en module 
à l'unité pour +, k,, k,, p, d quelconques. 

Ainsi le schéma aux différences (226), (227), tout comme 
le schéma implicite (217), est stable quelles que soient les 
relations existant entre les pas de réseau et, à la différence 
de (217), est économique du point de vue du nombre des 
opérations. En effet, le processus de résolution des équa- 
tions (226), (227) revient tout d'abord à la resolution du sys- 


tème (226) pour tout »7 donné, à la recherche de « et ensuite 
à la résolution du système (227) pour chaque # donné, ce qui 
donne w"*!, Les deux peuvent être résolus par la méthode 
habituelle du balayage. Les méthodes aux différences de ce 
type sont appelées méthodes des directions alternées ou 
méthodes des pas fractionnaires. 

Ilest facile decomprendre pourquoi la méthode du balayage 
matriciel se trouve être moins efficace que celle qui 
vient d’être exposée : elle est trop universelle. En effet, on 
sait que les avantages de la méthode habituelle du balayage 
sont dus à ce qu'on tient compte avec une grande précision 
de l'influence mutuelle des solutions dans les différents 
points. Considérons la méthode du balayage matriciel 
de ce point de vue. La relation entre les vecteurs uy 1, un 
écrite sous la forme (222) fait apparaître formellement les 
relations existant entre les composantes de ces vecteurs. 
ÏIl est cependant évident que l'influence mutuelle des diffé- 
rentes composantes diminue rapidement au fur et à mesure 
de l'éloignement réciproque des points de calcul corres- 
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pondants lors de l'augmentation de la différence entre les 
numéros 7. La méthode du balayage matriciel ne tient 
pas compte de cette particularité du système d'équations, elle 
vise une classe bien plus étendue de problèmes et dans le 
cas présent est loin d'être la meilleure. 

Nous avons envisagé seulement un seul des problèmes 
apparaissant lorsque l’on passe des problèmes unidimen- 
sionnels aux problèmes bidimensionnels. Tout comme pré- 
cédemment nous avons à cet effet utilisé un exemple parti- 
culicer caractéristique permettant de faire apparaître l'es- 
sence du problème. Il est évident qu'en augmentant le 
nombre des dimensions on voit apparaître de nouveaux 
problèmes, mais nous n’allons pas nous y arrêter. Ceux-ci 
sont essentiellement liés aux difficultés d'approximation des 
domaines multidimensionnels par des réseaux de calcul 
convenables et à l'obtention d'un algorithme de calcul 
simple et économique. 


Problèmes 


1. Trouver et étudier les différents schémas pour la 
solution de l'équation 


oU 1104 oU 

TT HR MAL PL a D 

en ulilisant la maille de calcul représentée fig. 16. 
Etudier également le schéma aux différences du typo 


un+î! n—i n n n n 
R, m UK, ue Unit. m Uk 1. m Lb point Uh, m1 = 
2t 2hx 2hy ès 


2. Etudier les méthodes itératives pour la résolution des 
systèmes d'équations (219), (220). Donner une estimation 
du nombre d'itérations. 

3. Donner une estimation du nombre d'opérations arith- 


méliques nécessaires à la résolution des problèmes sur un 
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intervalle de temps fini et dans un domaine limité, en 
utilisant 

a) le schéma explicite (214), 

b) le schéma implicite (217) résolu par la méthode du 
balayage matriciel, 

c) le schéma implicite (217) résolu par les méthodes 
itératives, 

d) la méthode des directions alternées (226), (227). 

Dans les trois derniers cas poser t = k. 

4. Eliminer entre les formules (226), (227) la grandeur 
intermédiaire uw. Comparer l'équation aux différences obte- 
nue avec le schéma implicite (217). 

5. Montrer que pour la solution du problème (213) on 
peut utiliser les algorithmes aux différences donnés par les 
formules 


Hé n 
Uh,m—Uk, m _ Uhai, me 2h, m + Un, m + 


T h2 
n n n 
them 28, m FER, m1 
Die ms à 
h} 
n+1i T7 n+i n+i 
Uhm = Uh, m __Ux, Va 1 20h on Um 
2 
T h? 
n n n 
Re m+1—2Ux, m'FUk, m—1 
h£ 
et 
La —u" Cond Crud nd 
Rom Uk, m _ Uhss, m— 2h, m+Uh-1, m 
T h2 L 
À. 7 +1 n+i{ n+i 
Uk tm UR, m x. ne 4 20 tm + Uh m 1 
T h£ 


Les comparer avec le schéma implicite (217) et la méthode 
des directions alternées (226), (227) (en excluant u). 
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6. Envisager la possibilité de généralisation de toutes 
les méthodes mentionnées dans le paragraphe présent aux 
problèmes tridimensionnels correspondants. 


$ 12. PROBLÈMES STATIONNAIRES 


On utilise ce terme pour les problèmes décrivant des 
états stationnaires des différents systèmes ne changeant pas 
dans le temps. À titre d'exemple typique on peut étudier le 
problème suivant. Soit à trouver la fonction U (x, y) satis- 
faisant dans un certain domaine limité G du plan x, y 
à l'équation 


D = (a 4), (229) 


prenant sur la limite l' de ce domaine des valeurs données, 
à savoir 


Ur =2£. (230) 


Il n'est pas difficile de trouver le problème aux diffé- 
rences correspondant. Recouvrons le domaine G par un 
réseau de calcul, pour plus de simplicité les pas suivant x 
et y sont pris égaux (fig. 20). Sur ce réseau on remplace 
l'équation (229) par la relation aux différences 


Un, m—2Uh, m+Uh-1, Un, mai — 2h, m4UR, m-1 
+ m or dd n + m+ ne = fh, is (231) 


ayant un sens pour chaque point interne de calcul. On en- 
tendra par point interne un point quelconque #, m pour lequel 
les quatre points voisins k + 1, m + 1 utilisés dans (231) 
sont disposés à l’intérieur du domaine G. Les autres points 


de calcul k, »m appartenant à G seront dits limites, leur 
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….. 


k,x 
Fig. 20 


ensemble étant désigné par y. On obtient les valeurs uy » 
sur y simplement par translation de la valeur g à partir du 
point voisin de la frontière 


u|,=£glr. (232) 


Le problème aux différences (231), (232) appartient à la 
classe envisagée au $ 5. Son étude revient à la vérification 
de l’approximation et de la stabilité. La première est évi- 
dente car à la limite, pour k — 0, (231) devient (229) et (232) 
devient (230) (la distance entre y et F étant de l'ordre de h). 

Nous allons démontrer la stabilité du problème aux diffé- 
rences établi, c’est-à-dire la coïncidence des ordres de gran- 
deurs de la solution et des seconds membres de (231), (232) 
quel que soit À. À cet effet on utilise l’artifice suivant. 
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Supposons que la solution &# du problème existe. Considé- 
rons les deux fonctions auxiliaires v, et v_ 


ve = Eu+a(r+ y) +h, (253) 


où «&,B sont pour le moment des constantes arbitraires. 
Désignons par Du le premier membre de (231). Substituons 
(233) dans (231), il vient 


Dv+ = +f + AG, 
car Du = f, D (x + y) = 4. DB — 0, 


On choisit & de telle sorte que partout dans le domaine G 
on ait l'inégalité 


Dv+ > 0. (234) 
Il est évident que pour cela il suffit de poser 
@= 7 max | fn. ml: (235) 
Puis on détermine B de telle sorte que 
VO. (236) 
En vertu de (232), (233) ceci aura lieu pour 
B — — Max | g|—œmax (a? + y?). (237) 


a 


Supposons que max v, corresponde à un certain point 
interne #,m. Dans ce cas au moins dans un des points 
voisins la valeur v, est inférieure à la valeur maximale, 
et comme on peut facilement le voir à partir de (231), 
(Du+}n, m < 0. Ceci se trouve en contradiction avec (234) 
et par conséquent max v, ne peut se trouver que sur la 
frontière y. Mais, vu (236), v+ y est négative, elle est donc 
négative partout, ainsi 


+ Ur, m + @ (2? + y}, m+B<0, 
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c'est-à-dire en vertu de (235), (237), on a 


max | uy, max [£ 147 max | fn m [max (4? + y). (238) 
k,m 

L'inégalité obtenue signifie que le problème (231). (232) 
est stable. 

Nous avons démontré également que la solution existe 
et qu'elle est unique. En effet comme chaque solution doit 
satisfaire à (238), alors pour g = f — 0 seule la solution 
triviale w — 0 est possible. Par conséquent la solution du 
système d'équations linéaires non homogènes (231), (232) 
existe et est unique. 

Abordons maintenant les méthodes de résolution du 
système d'équations (231), (232). Du point de vue historique 
les méthodes itératives sont apparues les premières. Nous 
allons décrire la plus simple d’entre elles. 

Résolvons chacune des équations (231) par rapport à la 
valeur u,, , au point central de la maille de calcul: 


Un, m =+ (UTEN m + U41, m + Up, m-1 me U, m+1 — hfy, m) (239) 


et utilisons cette formule pour les itérations. Le processus 
de calcul est très simple: sur chaque v-ième itération on 
calcule la moyenne arithmétique des valeurs 4%) ,41 aux 
points entourant le point central donné et l’on obtient 
l'approximation suivante uf”4). 

Etudions la convergence du processus. Posons 


(v) 
uÿ. m = Ur, m + 69 m 


OÙ Ux,m est la solution exacte du système (231), (232). 
Il est dans ce cas évident que l'erreur 6, ,, sera déterminée 
par le pee itératif suivant 


SALE 4 OQY m + +6, m + +1 + 5-1), 6tv+1) |, == 0 
(240) 
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Désignons par Ô le maximum du module de 6j”, nous 
allons raisonner comme suit. Comme 6f”,, est la moyenne 
arithmétique des quatre valeurs 621, m+1, | 4, [ne sur- 


passe pas ô, co qui est vrai pour tous les points sur toutes 
les itérations. Mais pour les points voisins des points limites 
on peut donner une estimation plus précise. À savoir, si ne 
serait-ce qu’un des points voisins de #, m#» est un point limite, 
où Ô — 0,en ce pointé, mon a 


3 — 
LUN 0458 — ô 


Cette dernière inégalité est vraie pour toute la couche 
limitrophe de points. Passons à la seconde itération, où 
l'influence de la frontière s'étend encore à une couche de 
points où 


og — = É5. 


Cette estimation est a fortiori vraie pour la première couche 
limitrophe de points. En continuant le raisonnement, on va 
se déplacer lors de chaque itération à l’intérieur du domai- 
ne G, pour les couches de points traversées on obtient alors 
l'estimation 


EURE < (1 ——)8. 


Enfin lors d’une certaine itération d'ordre n, n = 1/h, 
tous les points de calcul se trouvent épuisés. Ceci signifie 
que lors de x itérations l'erreur ô a diminué au moins de 
(1—4-") fois. Lors des n itérations suivantes, elle diminue 
encore du même nombre de fois, etc. Nous avons démontré 
que pour v —> co l'erreur Ê —+ 0, c’est à-dire que le pro- 
cessus itératif converge. 

Pour » itérations l'erreur diminue de (1—4-") fois, par 
conséquent pour une itération elle diminue en moyenne de 
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Fig. 21 


({1—4-")1/n fois, ou (comme x = 1/h) de (1—4-1/hjh © 
= À — h4-Uh fois. Généralement on caractérise la vitesse 
de convergence par la grandeur de la décroissance relative 
de l'erreur pour une itération, c’est-à-dire par le rapport 
(ô(V) — 6(v+1)) 6M — x. Dans le cas présent pour cette 
grandeur on obtient l'estimation suivante 


k— ha Ua, (241) 


En fait, dans de nombreux cas cette estimation est trop 
surhaussée, le processus itératif converge plus rapidement. 
Ainsi, si le domaine G est un rectangle (fig. 21), il est fa- 
cile d'obtenir une estimation plus précise. À cet effet remar- 
quons que la formule (240) décrivant l'évolution de l'erreur 
d'une itération à l'autre peut être interprétée comme une 
formule d'un problème aux différences 


g + 1) _ RS" 
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d'une classe que nous connaissons déjà ($ 6). La seule diffé- 
rence est que v est maintenant le numéro de l’itération et 
non le numéro de la couche temporelle. C'est pourquoi pour 
l'étude de l'évolution de l'erreur 6% on peut appliquer 
le critère spectral. Pour le moment on ne tient pas compte 
des conditions aux limites Ô], = 0 et l'on pose 


fn — Spectre + my), (242) 


Comme toujours, ô(+1 est égal à A6 et dans le cas pré- 
sent, on obtient facilement que 


Fee q + cos 1) 


+, (243) 


c'est-à-dire | À | & 1. Cette estimation nous assurait la sta- 
bilité des problèmes d'évolution, mais maintenant elle 
devient insuffisante. Seule l'inégalité rigoureuse | À | < Î 
nous convient, car ce n’est que dans ce cas que la conver- 
gence est garantie: Ô(" —+ O0 pour y —+ æ. Il est facile de 
voir que |A | = 1 correspond aux fonctions (242) qui ne 
satisfont pas aux conditions aux limites ô|, = 0 (celles-ci 
s'obtiennent pour @ = 1 = 0 ou @ = 1 = x), donc l’esti- 
mation |À | peut être précisée. 

Nous ne laisserons que les combinaisons des fonctions 
(242) qui s’annulent sur les limites de notre domaine rectan- 
gulaire G (fig. 21), c'est-à-dire satisfont à la condition 


Ôo, m—=ÔK, m=0, Ht= 0; 3:58 4e 
Ôp, 0 = Ôn, x =0, k= 0, 1, es 


Comme les exponentielles figurant dans (242) s'expriment 
en fonction de sin kp, cos Xp, sin m1, cos map, les fonctions 
nous intéressant sont des combinaisons de ces dernières. 
Pour satisfaire à la condition limite de gauche Ô,, , = 0, 
ces fonctions doivent contenir le facteur sin Aÿ. Si l’on 
exige que soit vérifiée la condition limite sur l'extrémité 
droite, pour À — X, on obtient l'égalité Xœ = 0. Ceci 
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n’est possible que pour Xœo = pr, où p est entier. Ainsi 
nous devons considérer seulement un ensemble discret 
et même fini 


Pr=P + p=1, 2, ..., K—1, (244) 


car pour les autres p on obtient les mêmes fonctions discrètes 
sin À», et pour p = 0 ou p = KX, une fonction nulle sur 
le réseau. 

Pour des raisons analogues nos fonctions ô,, doivent 
contenir le facteur sin #1, avec 


Va= 97 q=1, 2, ..., M—1. (245) 
Ainsi les fonctions discrètes 


Ôn, m = Sin Xp Sin Ma (246) 


pour PpŸ, quelconques définis par les égalités (244), (245), 
satisfont aux conditions aux limites. 

Substituons (246) dans (240) au lieu de SR Après des 
calculs simples on obtient : 


v+1) _ C08 Pp+C08 Wq 
Go RS 


m = sin Æ, sin MW, 


c'est-à-dire 64, »m (246) sont des fonctions propres de l'opé- 
rateur d'itérations, et les valeurs propres correspondantes 
sont données par la formule 

COS @n + COS 1, 


Ap, 4 = (247) 


coïncidant avec (243). La réserve des fonctions propres est 
grande (on peut montrer qu'elle est suffisamment grande) et 
la grandeur Àp, 4 (247) permet de juger de la vitesse réelle 
de la convergence du processus itératif. 

Les plus grandes valeurs de | À,, , | sont atteintes pour 
les valeurs limites ,,4,, c'est-à-dire pour |cos pp | = 
= cos (x/K) et | cos, | = cos (x/M1). Comme Xh et Afh 
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déterminent les dimensions du domaine G, c'est-à-dire qu'ils 
sont de l’ordre de l'unité, on peut écrire 


cos (URI cos GTA) = 


max | Ap,al= 5 7 


pq 
Par conséquent la caractéristique introduite ci-dessus de la 
vitesse de convergence des itérations est 


x — hi, (248) 


ce qui est évidemment mieux que (241) obtenu par une 
estimation grossière. 

Nous avons déjà noté l’analogie existant entre le pro- 
‘cessus itératif et le problème d'évolution aux différences. 
En fait, elle va bien plus loin. Chaque problème station- 
naire peut être considéré comme un cas particulier d'un 
problème d'évolution, où seul nous intéresse l'état final, 
stationnaire, et non le processus de l'établissement. Ceci 
influe également sur les méthodes de résolution des pro- 
blèmes stationnaires. Tous ces problèmes sont itératifs, 
à partir du plus simple: calcul des racines d'une équation, 
ils utilisent tous l’évolution, même fictive. Les problèmes 
linéaires font exception, mais si l'on tient compte du fait 
que même la division est une opération itérative, il devient 
évident que cette exception confirme la règle générale. Nous 
avons commencé par écrire l'équation f(x) — O0 sous la 
forme x—@ (x) ($ 1). Nous allons terminer en répétant cette 
procédure pour le problème du présent paragraphe. 

Ainsi pour la résolution des problèmes stationnaires on 
peut appliquer tout l'appareil élaboré pour les problèmes 
d'évolution. La formule envisagée du processus itératif (239), 
peut s'écrire sous la forme 


(v+1 \ (v 91N) (v) 
Um _ ui, … uv) 1, m —2u}, CR Uk}; .m 
h?;4 


h2 + 


TUE 1. 2), + TL 1 
+ D Re = În, m° 
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Si v est le numéro de la couche dans le temps, et h?/4, le 
pas t, on obtient le schéma aux différences (214) pour la 
résolution de l'équation bidimensionnelle de la conducti- 
bilité thermique que nous connaissons déjà. Remarquons 
que sa condition de stabilité se trouve vérifiée car t/x? — 1/4. 

Comme nous l’avons vu, la meilleure méthode de résolu- 
tion du problème mentionné s'est trouvée être la méthode 
des directions alternées. Cette méthode n'impose aucune 
restriction sur le pas + et il y a lieu de s'attendre à ce qu'en 
l'utilisant on atteindra plus rapidement l’état stationnaire 
limite, qui est la solution de notre problème (231), (232). 

C'est pourquoi nous allons revenir de nouveau aux for- 
mules (226), (227) en ajoutant le second membre /, h et 
nous allons les utiliser de paire avec les conditions aux 
limites correspondantes pour trouver la solution du problème 
stationnaire envisagé. Pour simplifier l’étude de la méthode, 
nous allons nous limiter au cas où le domaine G est un carré 
0 <Lzx, y L X et les pas du réseau sont égaux entre eux 
h = h, = = h La transition de la v-ième itération à la 
(v + {)-ième consiste à résoudre le système d'équations 


mr mn UN — Uni. m — Düy, m Hubs. m — a 
1/2 EE h2 
(v) 
u +u 
k.m+1— ufy m1 
T AUTRES nat, cn fn, m (249) 


4 


par rapport à uw, puis le système d'équations 


v+i Le. d Cd at 
u m ip, mi __ Uhsi, m—2Uh, mŸUh-1, m + 
T/2 h2 
el ),,(V+1) (v+1) 
u y —2u + u 
k, h, m h,m-—1 e 
de 5 — 2 — — fn,m (250) 


par rapport à utv+1), Les indices k, m prennent des valeurs 


de O à À — X/h, les grandeurs u, ut“+1) tout comme ul) 
sont données sur la frontière, c'est-à-dire pour k, m = 0, 
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K et à chaque point interne correspond une paire d'équa- 
tions (249), (250). 

Pour ce qui est de l'algorithme de résolution des systèmes 
(249) et (250), nous l'avons déjà étudié, c'est la vitesse de 
convergence du processus itératif qui nous intéresse mainte- 
nant. Cette étude est analogue à celle qui a été faite ci-dessus 
dans le cas d’une itération simple. Posant uW) = u + ô" 
on obtient pour l'erreur 8% le même problème aux diffé- 
rences (249), (250) avec f = 0 et des conditions aux limites 
nulles. Les fonctions discrètes 6,,,, (246) avec 


Pp = PIUK, Va = UK; p,q = 1,2, ..., K — 1 (251) 


seront de nouveau les fonctions propres de l'opérateur d’ité- 
rations. Les valeurs propres À déterminant la vitesse de 
convergence 8 “expriment par la formule (228) avec q = ,, 
Ÿ = Ÿ, c'est-à-dire elles sont des produits de facteurs de 
même type À = Anh, Ecrivons l'un d'eux 


Tu sin? — Tr 


Âp= SE Per (252) 
1 +7 sin DE 


Le second facteur À, s’oblient à partir du premier en rem- 
plaçant @» par Ÿ el par conséquent 


max|A[=max|A,|max|2,]— maxl]A,/. (253) 


Il est évident que max | À | < Î pour n'importe quelle 
valeur positive du paramètre 7, c'est-à-dire le processus ité- 
ratif converge toujours. Mais pour pouvoir estimer la vitesse 
de convergence, on aura besoin d'une estimation la plus 
exacte possible de la grandeur |A |. 

Comme Xh — X, p, a pour limites de variation 


hTUX Ep LA — hU/X, (254) 
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et pour À, pour À petit, on a l'estimation suivante (fig. 22) 


1—T— 


2)] ir —O (W2)] 
—— 5 ——— > — 5 —. (259) 
1+7T JE [1+ 0 (h°)] HT 7 [1— 0 (h?)] 


Nous avons besoin que max | À, | soit minimal. Nous 
avons à notre disposition le paramètre +, qui représente 
le temps fictif; il semble que plus grand est t, meilleure 
doit être la convergence, car on atteint alors plus rapide- 
ment l’état limite. En effet, lorsque t augmente, le premier 
membre de (255) décroît. Cependant, le second membre do 
(255) tendra alors vers —1, ce qui aura pour effet de ralentir 
la convergence. Cet effet peut s'expliquer par ce que, bien 
que pour t grand on se déplace plus rapidement dans le 
temps, ce déplacement est plus grossier, l'imprécision de la 
solution obtenue surpasse sa variation. 

Il est évident que la valeur optimale de t sera celle 
pour laquelle les premier et second membres de (225) sont 
égaux en module. En négligeant les infiniment petits 
d'ordre supérieur à © (k*), on obtient l'équation suivante 
pour T: 


n° 2 
tr 1-tir 
one 0 in dt 
HT 1+tr 


dont la solution est 


T = AX/n. (256) 


En substituant cette valeur de + dans l'inégalité (255), on 
obtient l'estimation 


max|A,|=1—Ahx/X +0 (h?), 
c'est-à-dire, en vertu de (253), 


max|A|—=1—2hn/X +0 (là), 
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Fig. 22 


et la grandeur x caractérisant la vitesse de convergence se 
trouve être de l'ordre de À 


2n 
x h—. (257) 

Ainsi le processus itératif utilisant la méthode des di- 
rections alternées avec T — À a une vitesse de convergence 
dix fois meilleure que la simple itération où x — h? (248). 

Si l'on développe l'erreur Ô en composantes qui sont 
des harmoniques de la fonction (246), (251), À LÀ, sera alors 
le coefficient d'amortissement de la composante pour une 
itération. La formule (252) et la fig. 22 montrent que les 
différentes composantes sont amorties différemment. Les 
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plus amortis sont les harmoniques de fréquences Çp pour 
lesquels 


Il est évident que le choix de t permet de régler la gammo 
des fréquences. Ainsi pour la valeur mentionnée de t (256) 
ce sont les fréquences moyennes. Mais leur amortissement 
trop important n'est pas intéressant car la vitesse de con- 
vergence est déterminée par les cocfficients d'amortisse- 
ment des fréquences extrêmes œ,, Va het pr, Ÿy an. 
Ceci porte à penser qu’il serait bon d'utiliser dans les 
différentes itérations différentes valeurs de +, ceci afin que 
d'assurer un amortissement uniforme de toutes les fré- 
quences. Ainsi, en utilisant des suites t — T+(% spécialement 
élaborées, on arrive à obtenir des méthodes dont la vitesse 
de convergence est encore plus grande que (257), par exemple 
(voir problème 2) 
1 

K — Mn 7%) : (258) 
Nous n'allons pas nous arrêter sur ces méthodes. 

Toutes les méthodes itératives de résolution des sys- 
tèmes d'équations aux différences (231) ont une vitesse de 
convergence faible qui décroît plus ou moins rapidement 
lorsque À augmente: x — 0 pour À — 0. Néanmoins elles 
sont plusavantageuses que lesméthodes directesnonitératives. 
Nous allons donner une estimation du nombre d'opérations 
arithmétiques Ÿ nécessaires à la résolution du problème. 

Commençons par la méthode générale d'exclusion. 
Lorsqu'on l'utilise, le nombre d'opérations nécessaires est 
de l'ordre du cube du nombre d'inconnues. Ce dernier étant 
de l’ordre de 1/k?, on a 

N = A/R, (259) 

Dans le paragraphe précédent nous avons décrit la 

méthode du balayage matricielle qui est, de toute évidence, 
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applicable dans le cas présent. Pour cette méthode il faut 
N = 1h (260) 


opérations arithméliques car pour inverser une matrice 
d'ordre {/X il faut 1/k% opérations, et on a en tout —1/h 
de telles matrices. 

Les méthodes itéralives donnent la solution approchée 
du système, c'est pourquoi il y a lieu d'estimer le nombre 
d'opérations nécessaires pour diminuer l'erreur d'un nombre 
donné de fois. Lors d’une itération l'erreur décroît de 1 — %x 
fois. Si l'erreur de l'approximation initiale est prise égale 
à l'unité, après x itéralions elle sera égale à 


(1—%x)"—=e. 
Par conséquent pour que la précision soit égale à e, il faut 
Ine In (1/e) 


D In ({ —%) é * 


itérations. Dans les méthodes envisagées le nombre d'opéra- 
tions arithmétiques pour une itération est de l'ordre du 
nombre de points du réseau, c'est-à-dire <1/h°. Par consé- 
quent, le nombre total d'opérations est 
In (1/e) 
di xh2 

En y substituant les différentes valeurs de x à partir 

de (248), (257), (258), on obtient : pour une itération simple 


N = In ({/e)/h*, (261) 
pour la méthode des directions alternées 
N = ln ({/e) h° (262) 
et dans le cas d’un choix spécial de + 
N = In ({/e) In (1/h)/h3, (263) 


En comparant (259), (260) avec (261) à (263), on voit 
l'avantage des premières. Il est évident que l’économie des 
méthodes itératives provient de ce que lorsqu'on les utilise 
on arrive à tenir compte au maximum des particularités du 
système d'équations. 
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Problèmes 


1. Déterminer la valeur optimale du paramètre v lors du 
processus itératif suivant les formules du type (249), (250) 
dans le cas où le domaine G est un rectangle de côtés X, 
Y et les pas À,, key ne sont pas égaux entre eux. 

2. Montrer qu'il est possible d'obtenir la vitesse de con- 
vergence des itérations donnée par la formule (258). À cet 
effet étudier le processus itératif (249), (250) à t variable: 


= zv-1, Vi: 2 san 
TO) = TV -n), v=n+ti, n+2, ..…. 


Choisir les valeurs 7% et » de telle sorte que la grandeur 


(v) 
Ë GE sin Es 


figurant dans l'expression de À, (252) pour tout p admissible 
se trouve pour un certain v — v (p) dans l'intervalle entre 


V z et 1/V/z. La valeur AM) correspondant à ce v doit vérifier 
l'inégalité 


(v} 1—V 3 

Le < 1+V3 

qui peut être utilisée pour estimer la valeur moyenne À, 

pour un cycle de » itérations et par conséquent la vitesse 
de convergence. 

3. Etudier la possibilité d'application des méthodes 
mentionnées dans le problème 5 du $ 11 au cas des problèmes 
stationnaires. 

4. Proposer un algorithme de calcul pour la solution de 
l'équation (229) dans les domaines de géométrie compliquée 
(figure formée de rectangles, cercle et autres) pour diffé- 
rentes conditions aux limites. 

5. Elaborer et étudier les schémas aux différences pour 
la résolution des systèmes d'équations 
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oU dv OV AU 
7 Ta —/@ y), Fr op END y) 


dans le rectangle 0 << 7 < À, 0 < y < Y avec les condi- 
tions aux limites sur les côtés 


U (0, y) = a (y), U(x, Y) =$ (a), 
V (x, 0) = y (a), V (4, y) = à (y). 


En particulier, envisager l'algorithme itératif utilisant 
des équations aux différences du type 


Canal Ph 


u — 
aus k,m-—1/à kR—1, m— 1/9 
Uk, m3; TT hr = 


: VR som UV, m1 
ns om EE 


ed Ph 
Cm dl 


: Vhaifa,m Vhetjam 
( R=1/2,m— EE 


hx 
uy, ma4l/ —u#, m—1/ 
=: Uh-1e, m—T a — TER. m: 
y 
vA nt, pris 
n+i 1/0. m “h-1/2, m1 
DR mg TE = 
k.m=—1/ h—3,m-1/2 
Un À Tfh 17e, m 1/9 
x 
un+i —urti 
n+1 h. m+1/ h,m-1/9 
DR m TE 2 
k+1/9,m h—1/9, m 
= Uh 1/9, m + nn TER, ms 
x 
k= 1, 2, 7 K; m = À, 2: .….., ÂT: 
X Y 
K + /2 M + 1/2 


où » est le numéro de l'itération, = désigne les valeurs 
intermédiaires des inconnues, + un certain paramètre. 
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